Compressive sampling and signal inference

ABSTRACT

A transmission mask or coded aperture is used in spectroscopy to compressively sample an optical signal. The locations of transmissive and opaque elements of the mask are determined by a transmission function. The optical signal transmitted by the mask is detected at each sensor of a plurality of sensors dispersed spatially with respect to the mask. A number of estimated optical signal values is calculated from sensor measurements and the transmission function. The optical signal is compressed by selecting the transmission function so that the number of measurements is less than the number of estimated optical signal values. A reconstructed optical signal is further calculated using signal inference. An imaging system created from plurality of encoded subimaging systems compressively samples an optical signal. Encoding methods include but are not limited to pixel shift coding, birefringent shift coding, transmission mask coding, micro-optic coding, diffractive coding, interferometric coding, and focal plane coding.

CROSS-REFERENCE TO RELATED APPLICATION

This application is a continuation application of U.S. patentapplication Ser. No. 11/183,838, filed Jul. 19, 2005, which claims thebenefit of U.S. Provisional Patent Application Ser. No. 60/589,475 filedJul. 20, 2004, U.S. Provisional Patent Application Ser. No. 60/589,530filed Jul. 20, 2004, and U.S. Provisional Patent Application Ser. No.60/622,181 filed Oct. 26, 2004. All of the above mentioned applicationsare incorporated by reference herein in their entireties.

BACKGROUND OF THE INVENTION

1. Field of the Invention

Embodiments of the present invention relate to systems and methods forcompressive sampling in spectroscopy and spectral imaging. Moreparticularly, embodiments of the present invention relate to systems andmethods for using less than one measurement per estimated signal valuein spectroscopy and spectral imaging.

2. Background Information

A signal is a physical phenomenon distributed over a space. Examplesinclude signals distributed over time, such as electromagnetic waves onantennas or transmission lines, signals distributed over Fourier space,such as optical or electrical spectra, and multidimensional signalsdistributed over physical space, such as 2D and 3D images. Digitalsignal analysis consists of signal estimation from discretemeasurements. For many years, sampling theory formed the theoreticalcore of signal analysis. Basic signal characterization required regularsampling at a rate proportional to the signal bandwidth. The minimalsampling rate is termed the Nyquist frequency.

Over the past couple of decades, however, novel sampling approaches haveemerged. One generalized sampling approach envisioned a bank of filteredcopies of a signal at sub-Nyquist rates. Since this generalized samplingapproach, a number of generalized sampling strategies and analyses haveemerged. Particular attention has been given sampling on multiresolutionbases and to irregular sampling. Multiresolution sampling is regarded asa generalization of generalized sampling and is used as a means ofbalancing sampling rates and parallel signal analysis.

Signal compression technology has developed in parallel with generalizedsampling theory. Multiresolution representations, implemented usingfractal and wavelet analysis, have been found to be critically enablingfor signal compression on the basis of the empirically observedself-similarity of natural signals on multiple scales. For appropriatebases, natural signals yield sparse multiscale representations. Sparsityand hierarchical self-similarity have been combined in signalcompression algorithms such as the embedded zero-tree and setpartitioning in hierarchical trees algorithms.

Generalized sampling and compression technologies resolve specificchallenges in signal processing systems. Generalized sampling enablessystems to sample at lower rates than the Nyquist frequency, and datacompression enables post-sampling reduction of the system data load.

Previous examples of sub-Nyquist sampling are divided into parameterestimation approaches and interpolation approaches. Parameter estimationuses sampled data to fit an a priori signal model. The signal modeltypically involves much greater constraints than conventional bandlimits. For example, one can assume that the signal consists of a singlepoint source or a source state described by just a few parameters. As anexample of parameter estimation, several studies have consideredestimation of the frequency or position of an unknown but narrow-bandsource from sub-Nyquist samples. Interpolation generates new samplesfrom measured data by curve fitting between known points.

In view of the foregoing, it can be appreciated that a substantial needexists for systems and methods that can advantageously provide forsub-Nyquist sampling that can be used in spectroscopy and spectralimaging.

BRIEF SUMMARY OF THE INVENTION

One embodiment of the present invention is a method for compressivelysampling an optical signal. An optical component with a plurality oftransmissive elements and a plurality of opaque elements is created. Thelocation of the plurality of transmissive elements and the plurality ofopaque elements is determined by a transmission function. The opticalcomponent includes but is not limited to a transmission mask or a codedaperture. The spectrum of the optical signal is dispersed across theoptical component. Signals transmitted by the plurality of transmissiveelements are detected in a single time step at each sensor of aplurality of sensors dispersed spatially with respect to the opticalcomponent. Each sensor of the plurality of sensors produces ameasurement resulting in a number of measurements for the single timestep. A number of estimated optical signal values is calculated from thenumber of measurements and the transmission function. The transmissionfunction is selected so that the number of measurements is less than thenumber of estimated optical signal values.

Another embodiment of the present invention is a spectrometer usingcompressive sampling. The spectrometer includes a transmission mask, adiffraction grating, a plurality of sensors, and a processor. Thetransmission mask has a plurality of transmissive elements and aplurality of opaque elements. The location of the plurality oftransmissive elements and the plurality of opaque elements is determinedby a transmission function. The diffraction grating disperses a spectrumof an optical signal across the transmission mask. The plurality ofsensors are dispersed spatially with respect to the transmission mask.Each sensor of the plurality of sensors detects signals transmitted bythe plurality of transmissive elements in a single time step. Eachsensor also produces a measurement resulting in a number of measurementsfor the single time step. The processor calculates a number of estimatedoptical signal values from the number of measurements and thetransmission function. The transmission function is selected so that thenumber of measurements is less than the number of estimated opticalsignal values.

Another embodiment of the present invention is an optical wavemeterusing compressive sampling. The optical wavemeter includes a slit, adiffraction grating, a mask, a complementary grating, and a lineardetector array. A monochromatic source is incident on the slit. Thediffraction grating produces an image of the slit in an image plane at ahorizontal position that is wavelength dependent. The mask has atwo-dimensional pattern of transmission variations. The mask producesdifferent vertical intensity channels for different spectral channels.The complementary grating produces a stationary image of the slitindependent of wavelength. The linear detector array measures verticalvariations in intensity of the stationary image. The mask is created sothat a number of measurements made by the linear detector is less than anumber of spectral channels sampled.

Another embodiment of the present invention is a method for temporallycompressively sampling a signal. A plurality of analog to digitalconverters is assembled to sample the signal. Each analog to digitalconverter of the plurality of analog to digital converters is configuredto sample the signal at a time step determined by a temporal samplingfunction. The signal is sampled over a period of time using theplurality of analog to digital converters. Each analog to digitalconverter of the plurality of analog to digital converters produces ameasurement resulting in a number of measurements for the period oftime. A number of estimated signal values is calculated from the numberof measurements and the temporal sampling function. The temporalsampling function is selected so that the number of measurements is lessthan the number of estimated signal values.

Another embodiment of the present invention is a method forcompressively sampling an optical signal using an imaging system. Theimaging system is created from a plurality of subimaging systems. Eachsubimaging system includes a subaperture and a plurality of sensors. Theoptical signal is collected at each subaperture of the plurality ofsubimaging systems at a single time step. The optical signal istransformed into a subimage at each subimaging system of the pluralityof subimaging systems. The subimage includes at least one measurementfrom the plurality of sensors of the subimaging system. An image of theoptical signal is calculated from the sampling function and eachsubimage, spatial location, pixel sampling function, and point spreadfunction of each subimaging system of the plurality of subimagingsystems. The sampling function is selected so that the number ofmeasurements from the plurality of subimages is less than the number ofestimated optical signal values that make up the image.

Another embodiment of the present invention is a method for spatiallyfocal plane coding an optical signal using an imaging system. Theimaging system is created from a plurality of subimaging systems. Eachsubimaging system includes an array of electronic detectors. Eachsubimaging system is dispersed spatially with respect to a source of theoptical signal according to a sampling function. The optical signal iscollected at each array of electronic detectors of the plurality ofsubimaging systems at a single time step. The optical signal istransformed into a subimage at each subimaging system. The subimageincludes at least one measurement from the array of electronic detectorsof the subimaging system. An image of the optical signal is calculatedfrom the sampling function and each subimage, spatial location, pixelsampling function, and point spread function of each subimaging systemof the plurality of subimaging systems. The sampling function isselected so that the number of measurements for the plurality ofsubimages is less than the number of estimated optical signal valuesthat make up the image.

Another embodiment of the present invention is a method for temporallyfocal plane coding an optical signal using an imaging system. Theimaging system is created from a plurality of subimaging systems. Eachsubimaging system includes an array of electronic detectors. The opticalsignal is collected at each array of electronic detectors of theplurality of subimaging systems at a different time step according to atemporal sampling function. The optical signal is transformed into asubimage at each subimaging system of the plurality of subimagingsystems. The subimage includes at least one measurement from an array ofelectronic detectors of each subimaging system. An image of the opticalsignal is calculated from the temporal sampling function and eachsubimage, time step, pixel sampling function, and point spread functionof each subimaging system of the plurality of subimaging systems. Thetemporal sampling function is selected so that the number ofmeasurements for the plurality of subimages is less than the number ofestimated optical signal values comprising the image.

Another embodiment of the present invention is a method for spatiallyand temporally focal plane coding an optical signal using an imagingsystem. The imaging system is created from a plurality of subimagingsystems. Each subimaging system includes an array of electronicdetectors. Each subimaging system is dispersed spatially with respect toa source of the optical signal according to a spatial sampling function.The optical signal is collected at each array of electronic detectors ofthe plurality of subimaging systems at a different time step accordingto a temporal sampling function. The optical signal is transformed intoa subimage at each subimaging system of the plurality of subimagingsystems. The subimage includes at least one measurement from an array ofelectronic detectors of the each subimaging system. An image of theoptical signal is calculated from the spatial sampling function, thetemporal sampling function and each subimage, spatial location, timestep, pixel sampling function, and point spread function of eachsubimaging system of the plurality of subimaging systems. The spatialsampling function and the temporal sampling function are selected sothat the number of measurements for the plurality of subimages is lessthan the number of estimated optical signal values that make up theimage.

Another embodiment of the present invention is a method forcompressively sampling an optical signal using multiplex sampling. Groupmeasurement of multiple optical signal components on single detectors isperformed such that less than one measurement per signal component istaken. Probabilistic inference is used to produce a decompressed opticalsignal from the group measurement on a digital computer.

Another embodiment of the present invention is method for compressivelysampling an optical signal. A reflective mask is created with aplurality of reflective elements and a plurality of non-reflectiveelements. The location of the plurality of reflective elements and theplurality of non-reflective elements is determined by a reflectionfunction. A spectrum of the optical signal is dispersed across thereflective mask. Signals reflected by the plurality of reflectiveelements are detected in a single time step at each sensor of aplurality of sensors dispersed spatially with respect to the reflectivemask. Each sensor of a plurality of sensors produces a measurementresulting in a number of measurements for the single time step. A numberof estimated optical signal values is calculated from the number ofmeasurements and the reflection function. The reflection function isselected so that the number of measurements is less than the number ofestimated optical signal values.

Another embodiment of the present invention is a method forcompressively sampling an optical signal using an optical component toencode multiplex measurements. A mapping is created from the opticalsignal to a detector array by spatial and/or spectral dispersion.Signals transmitted are detected by a plurality of transmissive elementsof the optical component at each sensor of a plurality of sensors of thedetector array dispersed spatially with respect to the opticalcomponent. Each sensor of the plurality of sensors produces ameasurement resulting in a number of measurements. A number of estimatedoptical signal values is calculated from the number of measurements anda transmission function, wherein the transmission function is selectedso that the number of measurements is less than the number of estimatedoptical signal values.

Another embodiment of the present invention is a spectrometer usingcompressive sampling. The spectrometer includes a plurality of opticalcomponents and a digital computer. The plurality of optical componentsmeasure multiple linear projections of a spectral signal. The digitalcomputer performs decompressive inference on the multiple linearprojections to produce a decompressed signal.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is an exemplary plot of a planar function space where a subspacespanned by the physical sampling is isomorphic to a subspace spanned bypossible signals, in accordance with an embodiment of the presentinvention.

FIG. 2 is an exemplary plot of a planar function space where a subspacespanned by the physical sampling is indicative of a subspace spanned bypossible signals, in accordance with an embodiment of the presentinvention.

FIG. 3 is a flowchart showing an exemplary method for compressivelysampling an optical signal using a spectrometer, in accordance with anembodiment of the present invention.

FIG. 4 is a schematic diagram of an exemplary spectrometer implementingarbitrary mapping, in accordance with an embodiment of the presentinvention.

FIG. 5 is a plot of an exemplary spectral transmission mask formultiscale Haar sampling, in accordance with an embodiment of thepresent invention.

FIG. 6 is an exemplary 8×8 transformation code for a quantized cosinetransform, in accordance with an embodiment of the present invention.

FIG. 7 is a comparison of exemplary simulated image reconstructionsusing permuted Hadamard transform, discrete cosine transform, andquantized cosine transform matrices using 4.69%, 15.63%, and 32.81% ofthe transformed components or available measurements, in accordance withan embodiment of the present invention.

FIG. 8 is a schematic diagram of an exemplary spectrometer employingcompressive sampling, in accordance with an embodiment of the presentinvention.

FIG. 9 is a plot of an exemplary spectral transmission mask used to makemeasurements on Haar-wavelet averages and details, in accordance with anembodiment of the present invention.

FIG. 10 is a flowchart showing an exemplary method for reconstructingHaar-wavelet averages and details for the exemplary spectraltransmission mask shown in FIG. 9, in accordance with an embodiment ofthe present invention.

FIG. 11 is a plot of an exemplary comparison of original neon spectraldata and reconstructed neon spectral data using the exemplary methodshown in FIG. 10, in accordance with an embodiment of the presentinvention.

FIG. 12 is a plot of an exemplary comparison of original neon waveletcoefficient data and reconstructed neon coefficient data using theexemplary method shown in FIG. 10, in accordance with an embodiment ofthe present invention.

FIG. 13 is a schematic diagram showing an exemplary optical wavemeter,in accordance with an embodiment of the present invention.

FIG. 14 is a plot of an exemplary mask for a 256-channel wavemeter withbinary sensors, in accordance with an embodiment of the presentinvention.

FIG. 15 is a plot of an exemplary mask for a 256-channel wavemeter with4-level sensors, in accordance with an embodiment of the presentinvention.

FIG. 16 is a depiction of an exemplary sample image taken using theexemplary optical wavemeter of FIG. 13 and the exemplary mask of FIG.14, in accordance with an embodiment of the present invention.

FIG. 17 is a plot of exemplary expected bit patterns for each of 100wavelengths using the exemplary optical wavemeter of FIG. 13 and theexemplary mask of FIG. 14, in accordance with an embodiment of thepresent invention.

FIG. 18 is a plot of exemplary measured bit patterns for each of 100wavelengths using the exemplary optical wavemeter of FIG. 13 and theexemplary mask of FIG. 14, in accordance with an embodiment of thepresent invention.

FIG. 19 is a plot of exemplary measured errors in bit patterns for eachof 100 wavelengths using the exemplary optical wavemeter of FIG. 13 andthe exemplary mask of FIG. 14, in accordance with an embodiment of thepresent invention.

FIG. 20 is a plot of exemplary quaternary digit patterns for each of 100wavelengths using the exemplary optical wavemeter of FIG. 13 and theexemplary mask of FIG. 15, in accordance with an embodiment of thepresent invention.

FIG. 21 is a plot of exemplary measured quaternary digit patterns foreach of 100 wavelengths using the exemplary optical wavemeter of FIG. 13and the exemplary mask of FIG. 15, in accordance with an embodiment ofthe present invention.

FIG. 22 is a plot of exemplary measured errors in quaternary digitpatterns for each of 100 wavelengths using the exemplary opticalwavemeter of FIG. 13 and the exemplary mask of FIG. 15, in accordancewith an embodiment of the present invention.

FIG. 23 is a flowchart showing an exemplary method for temporallycompressively sampling a signal, in accordance with an embodiment of thepresent invention.

FIG. 24 is a flowchart showing an exemplary method for spatially focalplane coding an optical signal using an imaging system, in accordancewith an embodiment of the present invention.

FIG. 25 is a flowchart showing an exemplary method for temporally focalplane coding an optical signal using an imaging system, in accordancewith an embodiment of the present invention.

FIG. 26 is a flowchart showing an exemplary method for spatially andtemporally focal plane coding an optical signal using an imaging system,in accordance with an embodiment of the present invention.

FIG. 27 is a flowchart showing an exemplary method for compressivelysampling an optical signal using multiplex sampling, in accordance withan embodiment of the present invention.

FIG. 28 is a flowchart showing an exemplary method for compressivelysampling an optical signal using a reflective mask, in accordance withan embodiment of the present invention.

FIG. 29 is a flowchart showing an exemplary method for compressivelysampling an optical signal using an optical component to encodemultiplex measurements, in accordance with an embodiment of the presentinvention.

Before one or more embodiments of the invention are described in detail,one skilled in the art will appreciate that the invention is not limitedin its application to the details of construction, the arrangements ofcomponents, and the arrangement of steps set forth in the followingdetailed description or illustrated in the drawings. The invention iscapable of other embodiments and of being practiced or being carried outin various ways. Also, it is to be understood that the phraseology andterminology used herein is for the purpose of description and should notbe regarded as limiting.

DETAILED DESCRIPTION OF THE INVENTION

The goal of compressive sampling is to reduce the number of measurementsneeded to characterize a signal. The defining characteristic ofcompressive sampling is that less than one measurement is registered inthe system per signal value estimated. Compressive sampling can beviewed as an implementation of data compression in the physical layer.

One embodiment of the present invention is a method for usingcompressive sampling in optical imaging and spectroscopy. Positivity isthe defining feature of images and spectra. The nominal detected signalis envisioned as a positive count of the number of photons detected overa certain spatial and temporal window. It is, of course, possible tointroduce negative and complex weights in optical measurements usingspatial, temporal or spectral filters. For optical imaging, compressivesampling enables cameras with more pixels in the estimated image than inthe electronic focal plane. Compressive sampling can be particularlyimportant in spectral and temporal sampling of image data, in whichcases the number of measurements per spatial resolution element can orcannot exceed one but the number of measurements per estimated sample inthe full space-time-spectral data cube is much less than one. Foroptical spectroscopy, compressive sampling enables spectralcharacterization with fewer detectors than conventional systems. Deeperspectroscopic significance arises, however, from using compressivesampling for direct molecular or atomic analysis without forming aconventional spectrum.

Linearity and predeterminism are the defining features of samplingstrategies. Linearity refers to the condition that measurements, g,consist of linear projections of a source state f. Formally themeasurement process is expressed as g=Hf+b, where b represents noise.Predeterminism refers to the condition that the measurementtransformation H is fixed prior to measurement and is independent of theoutcome of any single measurement. While there are many examples ofadaptive measurement in both imaging and spectroscopy, simplicitydictates focusing on single step non-adaptive measurement.

Compressive sampling integrates previous work in signal sampling,representation and compression, multiplex sensing and group testing.Compressive sampling is a novel approach for sub-Nyquist frequencymeasurement. One novelty is in the combination of generalized samplingand multiscale data compression strategies with two previously unrelatedapproaches: multiplex sensing and group testing. Multiplex sensingconsists of taking measurements that depend on multiple source points inthe object space under examination.

Multiplex sensing involves the concept of the object space. The objectspace terminology originates in spectroscopy, where Fourier and Hadamardtransform instruments are referred to as multiplex spectrometers incontrast to isomorphic instruments based on grating diffraction. Inoptical systems, multiplex sensing has also been applied in codedaperture and interferometric imaging systems. A variety of other sensormodalities are viewed as multiplex systems, including computedtomography and coherence tomography sensors.

Conventionally, the attractions of multiplex sensing have includedincreased throughput (e.g., better light detection efficiency) andimproved signal verses noise discrimination via weighing design. Intomographic systems, the additional advantage of native multidimensionalimage reconstruction arises. Recently reference structure tomography hasbeen explored as a new approach to multiplex sensing and it has beenshown that the additional advantage of compressive sampling can beobtained in multiplex systems.

Previously, it was shown that the number of measurements needed tocharacterize a sample can be reduced by adaptive multiscale sensing.This approach combined conventional measurements with prior knowledge ofthe expected multiscale structure of the object under analysis. Whilequite effective in guiding magnetic resonance image acquisition, thisapproach requires data be acquired in serial and that current and futuremeasurements depend on past measurements. These constraints areunattractive in real-time spectroscopy and imaging.

Compressive sampling extends the adaptive multiscale sensing approach byusing multiplex sensing to implement group testing on multiscalefeatures. Group testing has generally been applied to analysis ofdiscrete objects, such as coins, weights or chemical samples. Anotherembodiment of the present invention is a method that applies groupsampling to filter bank samples drawn from spatially distributedsignals. In this embodiment, group testing consists of making signalmeasurements which depend on disjoint groups of resolution or filterelements abstracted from the signal.

1. Sampling in Signal Representation and Measurement

Representation sampling is distinct from measurement sampling.Representation sampling refers to the mathematical representation of acontinuous function ƒ(x), x ε R^(d), in terms of finite or countablymany discrete numbers. A representation sampling system consists of arule for generating sample values from ƒ(x) and a rule for approximatingƒ(x) from the sample values. As an example, the sampling rule forNyquist-Shannon sampling returns values of ƒ(x) at equally spaced pointsand signal approximation consists of summing products of sample valueswith the sinc basis. Measurement sampling refers to the physical processof measuring “samples” of ƒ(x). For example, a camera measures samplesof the light intensity on a focal plane. The measurement samples are theraw signal values measured from the physical pixels.

Imagers and spectrometers infer representation samples from measurementsamples of physical signals. In many cases, as in a simple camera, it isassumed that representation samples are identical to measurementsamples. Of course, representation samples are not unique for a givensignal. Sensor systems are often modeled by continuous transformationsof ƒ(x), as in the Fourier or Radon transform, prior to discretesampling. In these cases, representations are constructed using discretebases on the transformed space, and measurement samples are treated asrepresentation samples on the transform space. This approach is taken,for example, in x-ray tomography. Representation samples can also takemany forms on the native continuous space of ƒ(x). In most cases it isassumed that this space is L²(R^(d)). Many bases spanning subspaces ofL²(R^(d)) have been considered over the past twenty years in the contextof wavelet and generalized sampling theories. Compression and multiscalesignal fidelity have been primary motivators in considering new bases,which have not generally been closely associated with measurementsampling.

Measurement sampling on transformed spaces are selected on the basis ofphysical characteristics of measurement systems. Tomographic systems,for example, measure Radon transforms rather than native samples becauseRadon data can be measured non-invasively. Generalized representationsampling, in contrast, has been considered primarily in post-processing.It is possible to design sensors for which measurements directlyestimate generalized representation samples, but to date examples ofthis approach are rare. In optical systems, physical constraints, inparticular, require that optical signals be non-negative, and modalcoherence constraints limit the capacity to implement arbitrary samplingschemes.

Isomorphism between representation samples and measurement samples isnot necessarily a design goal. In Radon and Fourier transform sensorsystems, the fact that measurement is not local on the native signalspace is generally treated as an annoyance. Generalized representationsampling focuses primarily on high fidelity description of the signal onthe native space. In contrast, the goal of compressive sampling is toensure that measured data is as useful as possible in estimating thesignal. Given a finite number of measurements, surprisingly littleanalysis of whether or not conventional samples represent a goodstrategy has been done. Issues of uniform verses nonuniform samplinghave been widely considered, but such analyses normally consider thesample itself to relatively simple and local in a continuous space.

A precise model for what is meant by a “measurement” and how themeasurement is associated with the signal is needed. The process ofmeasurement is associated with the process of analog to digitalconversion, and measurement is specifically defined to be the process ofabstracting a digital number from an analog process. A measurementsample is thus a digital number with at least partial dependence on ananalog state.

While representation samples are in theory continuous or discrete, inpractice representation sample values are also rational numbers. Theprimary difference between measurement samples and representationsamples is that representation samples are associated in some simple waywith signal values in a continuous space. Measurement samples are, inprinciple, associated with any property of the signal and need not belocally dependent on the signal value.

Optical sensing includes a long history of sensing based on non-localmeasurements, dating back to before the Michelson interferometer.Concepts of locality are also fluid in optics. For example, Fourierspace is the native signal space for spectroscopy. Over the past halfcentury, Fourier and Hadamard spectrometers have been particularlypopular non-local measurement systems. Sensors for which measurementsdepend on multiple points in the native signal space are calledmultiplex sensors in the optical literature.

Embodiments of the present invention examine and exploit deliberatelyencoded differences between representation samples and measurementsamples. In particular, these embodiments show that the number ofmeasurements needed to characterize a signal can be less than the numberof representation samples used to describe the signal. A system in whichless than one measurement is made per reconstructed representationsample is compressively sampled. Multiplex sensing is essential tocompressive sampling.

The process of transforming measurement samples to the representationspace is referred to as signal inference. In this section, basicrelationships between representation samples and measurement samples areestablished. The process of physical layer signal compression asimplemented by multiplex measurements is described and characterized.Linear and nonlinear decompression strategies are detailed. Finally,algorithms, designs, and performance of compressive sampling systems areconsidered in succeeding sections.

Representation sampling is described by sampling functions Ψ_(i)(x) andreconstruction functions φ_(i)(x) such that the signal f, f ε L²(R^(d)),is represented as $\begin{matrix}{{{f(x)} = {\sum\limits_{i}\quad{f_{i}{\phi_{i}(x)}}}},{where}} & (1) \\{f_{i} = {\left\langle {\psi_{i},f} \right\rangle = {\int{{f(x)}{\psi_{i}(x)}{{\mathbb{d}x}.}}}}} & (2)\end{matrix}$

In Whittaker-Shannon sampling, for example, φ_(i)(x)=sinc(π(x−i)) andΨ_(i)(x)=δ(x−i), assuming the sampling period is scaled to one. Arepresentation sampling model is a mapping between a sample space ofdiscrete vectors and a reconstruction space of continuous functions. Thefunction space associated with Shannon's sampling, for example, consistsof band-limited functions in L²(R), any function ƒ(x) in the subspacecan be exactly represented, or reconstructed perfectly from the samples,by Eqn. (1).

Discrete representation of a continuous function is the motivation forrepresentation sampling. Ideally, ƒ(x) lies within the subspace spannedby the sampling functions. In practice, the process of projecting ƒ(x)onto the sampling basis and representing on the reconstruction basis isa projection of ƒ(x) onto a subspace of the full range of possiblesignals.

For simplicity, the case that a model uses the same family oforthonormal functions φ_(i)(x) for both the sampling and thereconstruction functions is considered. Vector φdenotes the functionfamily and V_(c) denotes the associated function space. For any functionƒ(x) E V_(c), its orthogonal projection onto the space spanned byφ_(i)(x) is

φ_(i),f

φ_(i)(x).

φi, f

is referred to as the canonical component or representation off alongthe canonical axis φ_(i)(x). The canonical representation of f ε V_(c)is then the vectorf _(c)=∫_(x)(x)f(x)dx.  (3)

In general, f_(c) ^(T)φ is the projection of f into the space V_(c). Thecomponents of f_(c) are referred to as the canonical signal values. Inthe case of a digital imaging system, f_(c) is the digital image of theobject. In the case of a spectrometer, f_(c) represents the projectionof the spectrum of the source under test into the space V_(c).

In contrast to representation sampling, the goal of measurement samplingis to convert the analog realization of a signal into digital form. Themeasurement samples need not be simply connected to a reconstructionbasis and the measurement sampling functions need not be drawn from anorthogonal basis. A linear measurement system transforms ƒ(x) intodiscrete measurements, denoted by g, as followsg=∫ _(x) h(x)f(x)dx  (4)where h is the vector of measurement response functions h_(n)(x), whichis also referred to as the physical sampling functions. In contrast withthe sophisticated and complete basis of Shannon sampling, the rectangleaveraging function is common in measurement sampling. A digital camera,for example, integrates the value of the optical intensity over a squarepixel.

Measurement sampling is connected to representation sampling through thephysical sampling functions h_(n)(x). In designing an optical sensorsystem, h(x) and the representation basis φ_(n)(x) are selected suchthat h_(n)(x) ε V_(c). Then,h=Hψ  (5)withh _(ij) =∫h _(i)(x)ψ_(j)(x)dx.  (6)

In view of Eqn. (3) and Eqn. (5), Eqn. (4) can be expressed as followsg=Hf _(c) +b  (7)where b represents noise.

The goal of imaging systems and spectrometers is to estimate ƒ(x) viaf_(c) from the physical measurements g. Eqn. (7) is the basic model forthe forward transformation from ƒ(x) to measurements. Signal estimationconsists of inversion of this model. The utility and effectiveness ofthese systems is evaluated based on measures of the fidelity between theestimated and canonical signals. The measurement system is completelyspecified by matrix H. It is a multiplex system if H has non-zerooff-diagonal elements. A measurement system is compressive if the lengthof g is less than the length of f_(c) (i.e., the number of rows is lessthan the number of columns in H).

With compressive sampling, the null space of H is multi-dimensional.V_(h) represents the subspace of V_(c) spanned by the physical samplingfunctions. V_(⊥) represents the subspace orthogonal to V_(h) in V_(c).Functions in V_(⊥) produce zero measurement by the measurement system H.The canonical representation of function ƒε V_(c) can be decomposed intotwo components, f_(c)=f_(h)+f_(⊥), representing the components off inV_(h) and V_(⊥), respectively.

Since f_(⊥) produces zero measurement, Eqn. (7) reduces tog=Hf _(h) +b  (8)Eqn. (8) is generally well-conditioned for estimation of f_(h). Butestimation of f_(⊥) from g is impossible in general without furtherinformation.

Signal inference is the process of estimating f_(c) from f_(h) withadditional knowledge about the signals of interest. V_(f) denotes thesubset of V_(c) spanned by possible signals ƒ(x). If V_(h) is isomorphicto V_(f), g approximates f_(c). If V_(h) is indicative of V_(f), thereis a one-to-one mapping between g and f_(c).

The distinction between isomorphic and indicative measurement is shownin FIGS. 1 and 2, where V_(c) is a plane. FIG. 1 is an exemplary plot ofa planar function space 100, V_(c), where a subspace spanned by thephysical sampling 110, V_(h), is isomorphic to subspace spanned bypossible signals 120, V_(f), in accordance with an embodiment of thepresent invention. Isomorphic signal inference is attractive when V_(f)120 corresponds to a linear subset of V_(c) 100. In FIG. 1, V_(f) 120corresponds to a line on the plane. V_(h) 110 is made equal to line 120.

In FIG. 1, V_(f) 120 is a straight line, or an affine space in V_(c)100. A simple linear projection is be used to relate V_(h) 110 to V_(f)120. The projection completely describes the state of the signal f_(c).The sampling is therefore compressive and the inference is linear.

FIG. 2 is an exemplary plot of a planar function space 200, V_(c), wherea subspace spanned by the physical sampling 210, V_(h), is indicative ofa subspace spanned by possible signals 220, V_(f), in accordance with anembodiment of the present invention. Indicative signal inference isattractive when V_(f) 220 is indexed by a linear projection of V_(c)200. In FIG. 2, V_(f) 220 corresponds to a curve on the plane. If V_(h)210 is a line parallel to the vertical axis, a one-to-one correspondenceis drawn between points in V_(h) 210 and points in V_(f) 220. Each pointin V_(h) 210 is mapped on to the corresponding point in V_(f) 210 byadding a unique vector drawn from V_(⊥).

In FIG. 2, V_(f) 220 is a curve of one parameter. Its canonicalcomponents are dependent on the parameter. The component f_(h), bycompressive sampling, can be used as an index or as parameters of f_(c).When V_(f) 220 is known a priori, f_(h) uniquely determines the signalf_(c).

The challenge in isomorphic measurement is to discover the subset ofV_(c) that approximates V_(f). Discovery is accomplished byKarhunen-Loève decomposition, for example, and similar methods.Indicative measurement poses the dual challenges of characterizing V_(f)and of development of an algorithm to map g on to f_(c). Whileisomorphic inference is always linear in g, linear and nonlinear methodsof indicative inference can be considered.

The primary challenges in using signal inference to achieve compressivesampling are characterization of V_(f), design of H such that V_(h) isboth indicative of V_(f) and physically achievable in sensor hardware,and development of an efficient and accurate algorithm mapping g into anestimate for f_(h) and mapping the estimated f_(h) into an estimate forf_(c).

In the case of linear signal inference, V_(f) is typically characterizedby analyzing a representative sample of signals methods such asKarhunen-Loève decomposition. Estimation algorithms are relativelystraightforward in this case. In nonlinear inference, f_(h) is linearlydetermined by measurement g, but the relationship between g and f_(⊥) isnonlinear and represented via the integration of the information off_(h) with prior knowledge regarding V_(f). For example, efficient grouptesting sensors have been demonstrated for the case in which V_(f) isknown to consist of isolated point sources in V_(c). In an extreme case,it is known that f_(c) has at most one nonzero component. Signals ofmuch greater complexity for which an articulated description of V_(f) isunavailable are considered below.

For signal reference, various methods to explore and exploit therelationship between g and f_(⊥) are presented. For example, thecomponents of f_(c) from multi-scale group testing measurement have ahierarchically local relationship. After the estimation of f_(h),methods that attempt to optimize certain objective measures over thenull space of H are used, subject to certain constraints. Both theoptimization objectives and the constraints are based on the prioriknowledge of the targeted signals. The priori knowledge may be describedby physical models, statistical models, or their combinations. Forexample, one set of group testing methods rely on Bayesian inferencebased on statistical models of spectra and images as well as multiplexmeasurements.

2. Sampling in Spectroscopy and Imaging

Earlier it was shown that a linear sensor drawing measurement samplesfrom a continuous distribution is characterized by the representationsample to measurement sample mappingg=Hf _(c) +b  (9)

Here the physical realization of this mapping in optical spectrometersand imagers is considered.

In some optical systems, such as coherent holography, the optical fieldis characterized such that g, H and f_(c) are bipolar or complex. Inconventional photography and spectroscopy, however, f_(c) represents anonnegative power spectral density or radiance. It is assumed that thecomponents of g, H and f_(c) are nonnegative and real. The primarymotivation for this restriction is the realization that potentialbipolar or complex representations not withstanding, the vast majorityof images and spectra of interest are represented as nonnegative arrays.

Within the constraint of nonnegativity, the broad classes oftransformations H are implemented on spectra and images in opticalsystems. Many examples of Hadamard and Fourier transform spectroscopyand imaging systems have been implemented. Nonnegative Hadamard andFourier transforms are implemented in these systems by adding a constantpositive bias to the transformation kernel. Hadamard transformspectrometers, for example, sample the S matrix, where if H a Hadamardmatrix S=(1+H)/2. Similarly Fourier transform systems sample (1+R{F})/2where F is the Fourier transform kernel.

Selection of a physically achievable and numerically attractive kernel His the core of sensor system design. In some cases, such as 3D imaging,H must have a multiplex structure because identity mappings are notphysically achievable. In conventional spectroscopy and imaging,however, identity mapping is physically achievable if aperture size andsystem volume are not constrained. Despite the fact that the identitymapping is the only mapping with condition number one under thenon-negativity constraint, in some cases multiplexing is performed as ameans of achieving compressive sampling and of optimizing designconstraints such as sensor cost, sensor sampling rate, and designgeometry.

FIG. 3 is a flowchart showing an exemplary method 300 for compressivelysampling an optical signal, in accordance with an embodiment of thepresent invention.

In step 310 of method 300, an optical component with a plurality oftransmissive elements and a plurality of opaque elements is created. Thelocation of the plurality of transmissive elements and the plurality ofopaque elements is determined by a transmission function. The opticalcomponent includes but is not limited to a transmission mask or a codedaperture.

In step 320, the spectrum of the optical signal is dispersed across theoptical component.

In step 330, signals transmitted by the plurality of transmissiveelements are detected in a single time step at each sensor of aplurality of sensors dispersed spatially with respect to the opticalcomponent. Each sensor of the plurality of sensors produces ameasurement resulting in a number of measurements for the single timestep.

In step 340, a number of estimated optical signal values is calculatedfrom the number of measurements and the transmission function. Thetransmission function is selected so that the number of measurements isless than the number of estimated optical signal values.

In another embodiment of the present invention, the number of estimatedoptical signal values is calculated by multiplying the number ofmeasurements by a pseudo-inverse of the transmission function. Inanother embodiment of the present invention, the number of estimatedoptical signal values is calculated by a linear constrained optimizationprocedure. In another embodiment of the present invention, the number ofestimated optical signal values is calculated by a non-linearconstrained optimization procedure. In another embodiment of the presentinvention, the optical component is adaptively or dynamically variable.In another embodiment of the present invention, the adaptively ordynamically variable optical component is based on electro-opticeffects, liquid crystals, or micro mechanical devices. In anotherembodiment of the present invention, the optical component implementspositive and respective negative elements of an encoding.

Arbitrary nonnegative mappings are implemented in spectroscopy. FIG. 4is a schematic diagram of an exemplary spectrometer 400 implementingarbitrary mapping, in accordance with an embodiment of the presentinvention. Single mode spectrometer 400 is a baseline design. Inspectroscopy of a single spatial mode, the input field is mapped througha pinhole spatial filter 410. Pinhole 410 is collimated by sphericallens 420 and diffracted by grating 430. Grating 430 disperses thespectrum by creating a linear mapping between direction of propagationor wavenumber and wavelength. Cylindrical lens 440, following grating430, focuses the field along the dispersion direction such that spectralcomponents at coded mask 450 are separated into lines. If each line is acolumn, codes arranged perpendicular to the dispersion direction incoded mask 450 modulate the spectrum in rows. The ij^(th) pixel on codedmask 450 then corresponds to the ij^(th) component of H, h_(ij). Thevalue of h_(ij) equal to the transmittance of coded mask 450 at thecorresponding pixel and is set arbitrarily between zero and one. Thepower transmitted by each row is focused by cylindrical lens pair 460 onto single output detector 470, such that the set of measurements made onthe output detector is g=Hf, where f corresponds to the power spectraldensity of the input source. f is calculated from g and H usingprocessor 480, for example.

FIG. 4 demonstrates a form of compressive sampling in building anoptical wavemeter, in accordance with an embodiment of the presentinvention. While this approach demonstrates that any filter array can beimplemented in a spectrometer, the resulting system is physically nearlyas large as a fully sampled spectrometer with the same resolution.Another embodiment of the present invention is a multiplex spectrometerwith much smaller size that uses a photonic crystal filter. In thisembodiment, which is a compact system, it is not possible to implementarbitrary filter functions. In other embodiments of the presentinvention, coded mask 450 is implemented using a holographic filter, aholographic filter array, a multiplex grating hologram, a multipassinterferometer, a thin film filter, an array of thin film filters, andan array of photonic crystal filters.

Another embodiment of the present invention is a method involvingmultiscale sampling and group testing of correlated and uncorrelatedsignal components. Signals are random processes with finite correlationdistances. In the simplest model the signal component f_(i) is highlycorrelated with the signal component f_(i+1), but is independent off_(i+N) for sufficiently large N. Thus, if a group test is implementedon two uncorrelated coefficients linear combinations of f_(i) andf_(i+N) can be measured. Such delocalized group tests are easilyimplemented in a long baseline spectrometer system.

FIG. 5 is a plot of an exemplary spectral transmission mask 500 formultiscale Haar sampling, in accordance with an embodiment of thepresent invention. Transmission mask 500 is, for example, coded mask 450of spectrometer 400 in FIG. 4. In transmission mask 400, the spectrumconsists of 256 channels. The first 16 rows sample the mean spectrum oflength 16 blocks. The next 16 rows sample the means of the first eightchannels sampled in the previous blocks, the final 128 measurementssample individual channels. Mask 500 represents a unipolar nonnegativesampling of the Haar wavelet decomposition, with nonnegativity obtainedby zeroing the negative components as discussed above for Hadamard andFourier spectroscopy. Mask 500 is not compressive, 256 spectra arereconstructed from 256 measurements.

Arbitrary group testing is more challenging in imaging systems. Whileglobal transformations H are available in interferometric imagingsystems, it is much easier to implement mappings such that only acompact region of signal elements is active in any row of H,corresponding to measurements that sample only a local neighborhood inthe image. Pixel remapping for measurement in imaging systems isachieved through optical elements, such as multiple apertures,beamsplitters, gratings or birefringent elements or by electroniccombination of sensor values on the detection plane.

3 Linear Signal Inference

In linear signal inference, the estimated or reconstructed signal, viaits canonical representation f_(c), is linearly related to themeasurements g. Ideally, the sampling space V_(h) contains the sourcespace V_(s). While digital compression systems may implement essentiallyany transformation H for any specified source space, the range ofphysically achievable sensor transformations is more limited. In otherwords, V_(h) and V_(s) are not the same in many practical situations.The component f_(⊥) is estimated in the subspace orthogonal to V_(h),based on both the obtainable solutions and priori knowledge about thejoint space of V_(s) and V₁. The priori knowledge is described in twoaspects: the physical aspect and the statistical aspect. It can beexploited in signal estimate or reconstruction as well as in the designof sampling system H.

3.1 Compressive Sampling

The primary challenge of compressive sampling by linear inference liesin the joint optimization of the logical suitability and physicalplausibility of H. The design of sensing system H includes the selectionof representation basis functions Ψ and the design of samplingfunctions. In imaging, it is both desirable and feasible to selectrepresentational functions with spatially local support at multiplescales, such as a family of multi-scale wavelet functions. Inspectroscopy, a spectral function is expected to be continuous and haveone or more than one peak. Representational functions are selected sothat both the base and peak features are captured well with as few aspossible basis functions. Some representational functions, therefore,have global support.

3.1.1 Quantized Cosine Transform (QCT)

In one embodiment of the present invention, the elements of H are drawnfrom the set (−1,0,1). This embodiment is implemented usingcomplementary measurements H₁ and H₂ such that H=H₁−H₂ with elements ofH₁ and H₂ drawn from the binary set (0,1). Coding schemes based on suchmatrices can be implemented.

A non-compressive design is viewed as an extreme case where allmeasurements are used. In a compressive design, certain transforms areused to enable measurements of the principle components of the sourceimage in a representation and to enable source estimates by numericaldecompression with high fidelity.

In another embodiment of the present invention an image source ispartitioned into blocks of, for example, 8×8 pixels. Consider thedouble-sided transformation of each 8×8 block S_(ij), C_(ij)=QS_(ij)Q .FIG. 6 is an exemplary 8×8 transformation matrix 600, Q, for a QCT, inaccordance with an embodiment of the present invention.

Transform matrix 600 has the following properties. Its elements are fromthe set (0, 1, −1). The rows of matrix 600 are orthogonal. The rowvectors are nearly even in Euclidean length, with the ratio between thelargest and the smallest. When the source image is spatially continuouswithin block S_(ij), the transformed block C_(ij) exhibits thecompressible property that its elements decay along the diagonals. Theelements on the lower anti-diagonals are, therefore, truncated and theremaining elements are measured with fewer sensors. Denote by C _(ij)the truncated block matrix. An estimate of the source block S_(ij) isthen obtained from Q⁻¹ CijQ^(‡)(decompression). The same transformmatrix is used for all blocks of image S.

The above ideas are similar to the image compression with discretecosine transforms, as used in the JPEG protocol. In fact, the specificmatrix Q can be obtained by rounding the discrete cosine transform (DCT)of the second kind into the set (0, 1, −1). Q is referred to as thequantized cosine transform (QCT). But, the very structure of the QCTmatrix itself can be used to explain the compression.

FIG. 7 is a comparison 700 of exemplary simulated image reconstructionsusing permuted Hadamard transform (PHT), DCT, and QCT matrices using4.69%, 15.63%, and 32.81% of the transformed components or availablemeasurements, in accordance with an embodiment of the present invention.Visually the effectiveness of the compression with the QCT issurprisingly close to that with the DCT. A PHT matrix with row ordering[1, 5, 7, 3, 4, 8, 2, 6] is also used. Based on the basic 8×8 QCT andPHT matrices, larger transform matrices of hierarchical structure can beconstructed for multiple resolution analysis.

3.2 Signal Inference by Deterministic Characteristics

Another embodiment of the present invention is a method for linearsignal inference. More specifically, a first order approximation off_(⊥) from f_(h) is estimated based on priori knowledge of the sourcespace. Priori knowledge about deterministic features of a signaldistribution is described, accurately or approximately, by linearoperators mapping from f_(h) to f. The following optimization problemfor locating the deviation is based on the priori knowledge about thesource function $\begin{matrix}{{\min\limits_{f_{n} \in {{null}{(H)}}}{{\nabla\left( {f_{*} + f_{n}} \right)}}_{2}} = {\min\limits_{y}{{\nabla\left( {f_{*} + {Ny}} \right)}}_{2}}} & (10)\end{matrix}$where ∇ is the backward finite difference, or the discrete gradient, thecolumn of N span the null space of H, and ∥·∥₂ is the Euclidean vectornorm. In matrix form, the matrix for the backward finite difference,denoted by B_(∇), is a lower bidiagonal matrix with all the diagonalelements equal to 1 and all the supper-diagonal elements equal to −1.Via the normal equation corresponding to the minimization problem (10),the explicit linear dependence of f_(c) on f* is as followsf _(n) =Ny=N(C ^(T) C)⁻¹ C ^(T) B _(∇) f*  (11)where C=B_(∇)N. In practical computation, the matrix N can be obtainedonce for all from H such as via the QR factorization. The componentf_(n) is computed via a least-squares (LS) solution to (10) instead ofusing the closed-form expression of (11).

Another embodiment of the present invention is a method of determining alinear inference with minimal gradient. In step one, the input isobtained. The input includes sampling matrix H, measurement vector g. Instep two, a projection into the sampling space is calculated. Thiscalculation includes solving for a solution f* to the underdeterminedsystem of equations. In step three, a minimal gradient calculated. Thiscalculation includes obtaining the matrix N and solving f_(c) using theleast-squares solution in (10).

3.3 Signal Inference by Statistical Characteristics

Statistical knowledge of a signal distribution is exploited in at leasttwo ways. One is based on the Karhunen-Loève transform, the other isbased on conditional probabilities. The statistics are drawn from asufficiently large data set. The data set is represented with a matrixX. The j-th column corresponds to the samples of the j-th randomvariable x_(j). S is the matrix obtained from X with the mean translatedto zero and the variances normalized to 1.

3.3.1 Inference Based On The KL Transform

The discrete KL-transform, also called the Hotelling transform, isdetermined as follows. S=UΣV^(T) is the singular value decomposition ofS, where U and V are orthonormal and Σ is diagonal and nonnegative. Thediagonal elements, σ_(j)=Σ(j,j), are referred to as the singular valuesof S and they are nonincreasingly ordered. The transformed variables inthe transformed data Y=SV are uncorrelated and ordered withnonincreasing variances σ_(j) ². This is verified by checking thatY^(T)Y=Σ² is diagonal. It is also possible to verify that thetransformed variables remain zero means. The KL-transform is appliedoften for the so-called Principal Component Analysis (PCA) where thenumber of decorrelated variables with dominant variances (i.e., theprincipal components Y_(p)) is significantly smaller than the number ofthe original variables.

The samples of a source distribution f, with the mean shifted to zero,is represented more compressly in the principal components, Y_(p)f_(p).H is designed to determine the principal components. The principalcoefficients f_(p) are then determined by the following equationg=(HY _(p))f _(p).  (12)

The system of equations in (12) is no longer underdetermined and is wellconditioned, (i.e., f_(p) can be uniquely determined and stablycomputed). To get the mean value of the samples, f_(m), H is designed sothat the total sum of the samples can be measured, which is feasible inmany practical situations. The samples are then estimated byY_(p)f_(p)+f_(m).

4. Nonlinear Signal Inference

Nonlinear signal inference consists of estimating f_(⊥) by combining themeasured value of f_(∥) with prior knowledge of the signal subspaceV_(f). Formally, the estimated signal isf _(est) =f _(∥est) +f _(⊥est)(f _(∥est))  (13)As discussed earlier, nonlinear inference is attractive when thereexists a one-to-one mapping between f_(c) and f_(∥) even though V_(f) isnot a linear projection of V_(c). In this case, f_(∥) is completelyindicative of f_(c) and f_(⊥est) (f_(∥est)) is the point in V_(⊥)indicated as a component of f_(c).

The primary challenges of nonlinear inference are characterization ofV_(f), discovery of a linear mapping H such that V_(∥), is indicative ofV_(f), and development of an algorithm associating each point in V_(∥),with a corresponding point in V_(f). While solutions to these challengeshave been demonstrated in highly constrained situations, such aswavemeters and tracking systems, exact solutions are unavailable ingeneral spectroscopy and imaging systems.

Specification of f_(⊥est) (f_(∥est)) is challenging because for broadnatural signal classes, such as optical spectral and images, it is notpossible to precisely describe V_(f). One embodiment of the presentinvention is a method of overcoming this challenge using Bayesianinference. Bayesian inference is based on statistical characterizationof V_(f). Distributions for components of f given measurement states areestimated from a library of known signals. These distributions are thenused to determine likely signal states given measurements from unknownsignals. Bayesian inference produces an estimate corresponding to themost likely value of f_(c) given the measurement g. System design forBayesian inference consists of choosing H such that the most likelysignal corresponds to the actual signal for a given measurement. Inpractice, global maximization of the probability given a measurementstate is less practical.

Another embodiment of the present invention is a method using objectiveinference. Objective inference is based on plausibility measures off_(est). Plausibility is expressed in one or more “objective functions,”θ(f_(est)). Objective functions are measures of signal continuity,smoothness, or coherence. Objective inference selects f_(⊥est)(f_(∥est)) to maximize objective measures of f_(est). The challenge ofobjective inference is to discover a mapping H and an objective functionsuch that f_(est) is a good estimate of f_(c).

Another embodiment of the present invention is a method using graphicalinference. Graphical inference constructs objective functions fromempirical samples drawn from V_(f). By constructing empiricaldistributions and objective measures graphical inference is viewed as anintegration of Bayesian and objective methods.

5. Compressive Spectroscopy

Physical implementations of compressive sampling include but are notlimited to spectroscopy, digital imaging and temporal signal analysis. Aspectrometer measures a set of “spectral projections” of the formm[n]=S(λ),h_(n)(λ)

.  (14)Multiplex spectroscopy under which h_(n) corresponds to a Fouriercomponent or a Hadamard/Walsh code is common. More generally, h_(n)represents a filter function. Almost any real positive valued filterfunction can be implemented. Bipolar and complex filters are constructedfrom linear transformations of real positive filters.

FIG. 8 is a schematic diagram of an exemplary spectrometer 800 employingcompressive sampling, in accordance with an embodiment of the presentinvention. Spectrometer 800 includes spectral source 810, diffractiongrating 820, lens 830, mirror 840, mask 850, complementary mirror 860,complementary lens 870, complementary diffraction grating 880, anddetector array 890. Unknown spectral source 810 is imaged onto a slit atthe input. Spectral source 810 is assumed to uniformly illuminate theslit. If spectral source 810 is a point source, a combination ofcylindrical lenses can be used to expand spectral source 810 to fullyilluminate the slit normal to the plane of FIG. 8 while keeping spectralsource 810 focused in the plane of the FIG. 8. The input slit is imagedonto mask 850 through diffraction grating 820 or a hologram. The effectof diffraction grating 820 is to separate the image by color component.In spectrometer 800, blue is imaged at the top of the mask, green in themiddle, and red at the bottom.

Mask 850 is divided into rows. Each element in a row is a filtercomponent, modulating the color channel corresponding to thecorresponding column. The modulation is by absorption, meaning that thefilter weights are all positive. Negative or complex filter weights areobtained by linear weighting of the output filtered signals. When theslit is re-imaged through complementary grating 880 onto output detectorarray 890, each element on the detector array consists of a filteredprojection of the input spectrum. Since the filter function correspondsexactly to the transmittance of the mask, any filtering operation isimplemented by spectrometer 800. Thus, it is possible to implement acompressive sampling scheme in spectrometer 800.

5.1 Haar-Wavelet Based Spectral Sampling

Another embodiment of the present invention is a method for compressivesampling using Haar-wavelet based spectral sampling. The target spectraconsists of discrete projections $\begin{matrix}{{s\lbrack n\rbrack} = \left\langle {{S(\lambda)},{{rect}\left( \frac{\lambda - {n\quad{\delta\lambda}}}{\delta\quad\lambda} \right)}} \right\rangle} & (15)\end{matrix}$

The resolution of spectrometer 800 is δλ and there are N spectralchannels. So the slits in mask 850 are designed to pass spectral signalsof width δλ on to detector array 890. The slits are of discrete width sothat measurements on spectrometer 800 take the form $\begin{matrix}{{\overset{\_}{m}\left\lbrack n^{\prime} \right\rbrack} = {\sum\limits_{n}\quad{\alpha_{n^{\prime}n}{s\lbrack n\rbrack}}}} & (16)\end{matrix}$

Using a binary masks, σ_(n′n) assumes values 0 or 1. However, withineach row (corresponding to the n_(th) measurement) any sequence σ_(n′n)can be implemented. In another embodiment of the present invention andin order to maximize throughput, Hadamard S-Matrix codes of s[n′l] areimplemented in each mask row. For example, $\begin{matrix}{{s\left\lbrack {n^{\prime},l} \right\rbrack} = {\sum\limits_{n = {{{({n^{\prime} - 1})}l} + 1}}^{n = {n^{\prime}l}}\quad{s\lbrack n\rbrack}}} & (17)\end{matrix}$

Compressively sampled spectroscopy consists of making measurements ofthe form $\begin{matrix}{{m\left\lbrack n^{\prime} \right\rbrack} = {{\sum\limits_{n = 1}^{l{({J,N})}}\quad{{a_{J}\left\lbrack {n^{\prime},n} \right\rbrack}{{fa}_{J}\lbrack n\rbrack}}} + {\sum\limits_{j = 1}^{J}\quad{\sum\limits_{n = 1}^{l{({j,N})}}\quad{{d_{j}\left\lbrack {n^{\prime},n} \right\rbrack}{{fd}_{j}\lbrack n\rbrack}}}}}} & (18)\end{matrix}$where fa_(J)[n] describes the signal averages at level J using aHaar-wavelet transform, fd_(J)[n], 1≦j≦J describes the Haar details forlevels 1 through J. In another embodiment of the present invention, thefollowing measurements are made. At level J, 2 ^(−J)N averages aremeasured (i.e., a_(J)[n′n]=δ_(nn′)). At level J, 2 ^(−J)N details arealso measured (i.e., d_(J[n′,n]=δ) _(nn′)). For levels 1≦j≦J, 2^(−J)Nvalues of sums of some permutations of details are measured, (i.e.,d_(j)[n₁n]=Σ_(n) _(p) ₌₁ ^(2(j−j)δ) _(nf(n,n) _(p)) where f(n,n_(p)) isa perturbation of signal positions n′>N).

FIG. 9 is a plot of an exemplary spectral transmission mask 900 used tomake measurements on Haar-wavelet averages and details, in accordancewith an embodiment of the present invention. Mask 900 measures 256spectral channels with just 80 measurements, implementing 31.25%compression.

Mask 900 makes measurements using a “lifting” procedure. A Haar-wavelettransform is implemented with the following direct transform:$\begin{matrix}{{A_{J + 1}\lbrack n\rbrack} = \frac{{A_{J}\left\lbrack {{2n} - 1} \right\rbrack} + {A_{J}\lbrack n\rbrack}}{\sqrt{2}}} & (19) \\{{D_{J + 1}\lbrack n\rbrack} = \frac{{A_{J}\left\lbrack {{2n} - 1} \right\rbrack} - {A_{J}\lbrack n\rbrack}}{\sqrt{2}}} & (20)\end{matrix}$where A_(j) and D_(j) are the average and detail coefficients at levelj. The inverse transform is given by $\begin{matrix}{{A_{J}\left\lbrack {{2n} - 1} \right\rbrack} = \frac{{A_{J + 1}\lbrack n\rbrack} + {D_{J + 1}\lbrack n\rbrack}}{\sqrt{2}}} & (21) \\{{A_{J}\left\lbrack {2n} \right\rbrack} = \frac{{A_{J + 1}\lbrack n\rbrack} - {D_{J + 1}\lbrack n\rbrack}}{\sqrt{2}}} & (22)\end{matrix}$

Mask 900 measures all 32 average coefficients at level 3. The Haaraverage coefficients are the sum of spectral channels of appropriatelength and scaled to an appropriate factor depending on the level. Theseaverages at level 3 are used to find the averages and details at level 4using the direct transform defined in Eq. (19).

For the remaining measurements, measurements of sums of combinations ofdetails for levels 1 through 3 are made. Specifically, sums of 2 level 3details, sums of 4 level 2 details and sums of 8 level 1 details aremeasured. Since averages can be easily measured with mask 900, sums ofthe averages at each level are measure and the sums of details are thenestimated with the prior knowledge of all the parent details and parentaverages.

Consider the measurement m_(j)[n] produced by sensor n at level j givenby $\begin{matrix}{{m_{j}\lbrack n\rbrack} = {\alpha_{j}{\sum\limits_{k}\quad{A_{j}\lbrack k\rbrack}}}} & (23)\end{matrix}$where k includes all the averages selected by the mask for themeasurement at n. From the definition of the Haar transform, theproportional factor α_(j) at level j is given byα_(j)=√{square root over (2)}^(j)  (24)The averages at level j are related to the averages and details at Levelj+1 as given in Eq. 21. Therefore, $\begin{matrix}{\frac{m_{j}\lbrack n\rbrack}{\alpha_{j}} = {{\sum\limits_{k^{\prime}}\quad{A_{j + 1}\left\lbrack k^{\prime} \right\rbrack}} + {\sum\limits_{k^{\prime}}\quad{D_{j + 1}\left\lbrack k^{\prime} \right\rbrack}}}} & (25)\end{matrix}$where k′ represents the parent averages and details corresponding toindex k. $\begin{matrix}{{\sum\limits_{k^{\prime}}\quad{D_{j + 1}\left\lbrack k^{\prime} \right\rbrack}} = {\frac{m_{j}\lbrack n\rbrack}{\alpha_{j}} - {\sum\limits_{k^{\prime}}\quad{A_{j + 1}\left\lbrack k^{\prime} \right\rbrack}}}} & (26)\end{matrix}$

The sum of the details in Eq. 26 is then assigned to each of theindividual detail coefficient based on a weighting scheme. Specifically,a proportional weighting scheme is used in all the reconstructions.

FIG. 10 is a flowchart showing an exemplary method 1000 forreconstructing Haar-wavelet averages and details for mask 900 shown inFIG. 9, in accordance with an embodiment of the present invention.Method 1000 can be used for any mask.

FIG. 11 is a plot of an exemplary comparison 1100 of original neonspectral data and reconstructed neon spectral data using the exemplarymethod 1000 shown in FIG. 10, in accordance with an embodiment of thepresent invention.

FIG. 12 is a plot of an exemplary comparison 1200 of original neonwavelet coefficient data and reconstructed neon coefficient data usingthe exemplary method 1000 shown in FIG. 10, in accordance with anembodiment of the present invention.

5.1 Measurement Efficient Optical Wavemeters

Another embodiment of the present invention is a multiplex opticalsensor that efficiently determines the wavelength of aquasi-monochromatic source. The measurement efficient optical wavemeteris a special case of a general spectrometer based on compressivesampling in that the source is assumed to have only one spectralchannel.

5.2.1 Theory

A digital optical spectrometer estimates features of the spectraldensity, S(λ) of a source from discrete measurements. Measurementsconsist of spectral projections of the formm[n]=∫h _(n)(λ)S(λ)dλ,  (27)where n is the index of the discrete array of measurements. Measurementsare taken at different points in time (e.g. serially) or at differentpoints in space (in parallel), or in some spatial and temporalcombination. Measurements taken at different points in time rely on oneor just a few optical detectors, with changes in the optical systemaccounting for variation in h_(n)(x) from one measurement to the next.Measurements taken in parallel require an array of optical detectors.

An optical wavemeter characterizes S(λ) under a prior constraint thatonly one spectral channel is active. This assumption is expressed asS(λ)=αδ(λ−λ_(o)),  (28)where α is the amplitude of the signal at λ=λ_(o). In this case, themeasurements are:m[n]=αh _(n)(λ_(o))  (29)

The goal of the wavemeter is to estimate α and λ_(o) from themeasurements m[n]. The number of possible outcomes from the estimationprocess is D_(α)N_(λ), where D_(α) is the dynamic range estimated in α,and N_(λ) is the number of possible values of λ_(o). The number ofpossible outcomes from the measurement process is D_(m) ^(N) ^(m) ,where D_(m) is the measurement of dynamic range, and N_(m) is the numberof measurements. If the sensor system is designed to achieve aone-to-one mapping between measurement outcomes and sensor outcomes, thenumber of measurements required is $\begin{matrix}{N_{m} \geq \frac{{\log\quad D_{\alpha}} + {\log\quad N_{\lambda}}}{\log\quad D_{m}}} & (30)\end{matrix}$

As an example of a maximally efficient sensor, suppose that D_(α)=D_(m).For a monochromatic source, α is measured using a single measurementwith h₁(λ)=1. log N_(λ)/log D_(m) additional measurements are needed todetermine λ_(o). If N_(λ)≦log D_(m) then only one additional measurementis needed. If h₂(λ)=K(λ−λ_(min)) is selected for the second measurement,where K is a constant, then:m[2]=ακ(λ−λ_(min))  (31)

α is already known from m[1] and the dynamic range is assumed sufficientto descriminate N_(λ) spectral bins in the range from m[2]=0 tom[2]=ακ(λ_(max)−λ_(min)).

While the method described by Eq. (31) shows that it is possible tobuild an optical wavemeter with just two measurements, it alsoillustrates that high dynamic range and a carefully designed andcharacterized sensor spectral response are necessary to achieve suchefficiency. In practice, it is difficult to implement a filter-detectorcombination with perfectly linear spectral response. Efficient wavemeterdesign requires specification of both the dynamic range of themeasurement system and control over the spectral responses of themeasurements.

Suppose that h_(n)(λ) is independently specified within each wavelengthchannel to a dynamic range of D_(h). Possible values for h_(n)(λ) rangefrom 0 to h_(n max). Assuming that the dynamic range of the measurementsystem is sufficient for all values of h_(n), Eq. (29) produces D_(h)^(N) ^(m) different measurement states. Assuming that a has beencharacterized by a single measurement as before, N_(m−1) additionalmeasurements characterize λ₀ if D_(h) ^(N) ^(m−1) ≧N_(λ), orN _(m)≧1+log_(D) _(h) N _(λ)  (32)

Another embodiment of the present invention is a method for achievingthe limit specified in Eq. (32). The method assumes that the N_(λ)possible outcomes for estimation of λ_(o) are ordered and assignedindices n(λ) such that 1≦n(λ)≦N_(λ). This ordering is such that n(λ) ismonotonic in λ−λ_(min), but such proportionality is not necessary andphysical design constraints may make some other ordering attractive. Letd(n, l, b) be the l^(th) digit of n in the base b. A filter array suchthat $\begin{matrix}{{h_{n}(\lambda)} = {{d\left( {{n(\lambda)},n,D_{h}} \right)}\frac{h_{n\quad\max}}{D_{h}}}} & (33)\end{matrix}$enables characterization of λ_(o) in N_(λ)/log D_(h) measurements. Themethod has been described in Eq. (31) for the case D_(h)=N_(λ), yieldingN_(m)=2. In the more general case, comparing of Eqs. (29) and (33) itcan be seen that m[n]/(αh_(n max)) is the n^(th) digit in the base D_(h)representation of n(λ_(o)).

Under the method described earlier, h_(n)(λ) is programmed usingtransmission masks in a grating-based Fourier filter system. Thetransmittance of such masks is normally determined in a multiple stepcoating process. The number of transmittance values is typicallyD_(h)=2^(steps), where steps is the number of deposition steps. With 5deposition steps, for example, Eq. (32) suggests that a 1024 channelsource can be characterized using just 3 measurements. With just onedeposition step, a 1024 channel source is characterized in 11measurements.

Exemplary demonstrations of optical wavemeters for D_(h)=2 and D_(h)=4are provided.

5.2.2 Design

As described above, a key step in the compressive wavemeter is theimprinting of a spatial-spectral intensity code on the source. Asubsequent measurement of the spatial variation of the intensity,coupled with knowledge of the spatial-spectral code, then allowsreliable estimation of the input spectrum.

FIG. 13 is a schematic diagram showing an exemplary optical wavemeter1300, in accordance with an embodiment of the present invention.Wavemeter 1002 includes source 1305, slit 1310, lens 1315, diffractiongrating 1320, lens 1325, intensity encoding mask 1330, lens 1335,diffraction grating 1340, lens 1345, and linear detector array 1350.Wavemeter 1300 utilizes a one-dimensional spatial-spectral intensitycode. In other words, 1-D positional variation in intensity is producedthat is a function of wavelength.

This encoding is not produced directly, but instead as the result of athree stage process. First, diffraction grating 1320 converts thespectral density of source 1305 into a spatially-varying intensityprofile. Because source 1305 is monochromatic and grating 1320 islocated in the center of a 4-f imaging system, this produces an image ofthe input slit in the image plane at a horizontal position that iswavelength dependent.

Second, mask 1330 with a 2-D pattern of transmission variations isplaced in the image plane. The image of slit 1310 interacts with anarrow horizontal region of mask 1330 and is intensity-modulated by thevertical transmission pattern of mask 1330 at that horizontal location.Because the horizontal image location is wavelength dependent, ahorizontal variation in the mask transmission produces differentvertical intensity patterns for different spectral channels.

Finally, diffraction grating 1340 and imaging setup undo the correlationbetween spectral channel and horizontal position, producing a stationaryimage of slit 1310 regardless of wavelength. The vertical variation ofintensity that was imprinted by mask 1330, however, remains. Byarranging linear detector array 1350 along the vertical direction of theslit image, this intensity variation can be measured. Using theknowledge of the transmission pattern on mask 1330 and the dispersionproperties of grating 1320 and grating 1340, the spectral channel ofsource 1305 is uniquely determined.

In another embodiment of the present invention, linear detector array1350 is a two-dimensional array or a single element detector. In anotherembodiment of the present invention, the linear detector array 1350 is abase-N encoded mask and measures N levels of intensity. In anotherembodiment of the present invention, mask 1330 is a single elementphotodiode and is dynamically variable.

Another embodiment of the present invention is a method of creatingtransmission mask 1330. As described earlier, if it is assumed thatdetectors 1350 have m distinguishable output levels, then the optimalcode for N spectral channels requires 1+log_(m)N measurements. Thesemeasurements determine not only the spectral channel of the source, butalso its power. (For the special case of m=2, the spectral channel alonecan be determined by log₂N measurements, but a similar reduction is notpossible for other values of m.)

If N spectral channels are numbered from 1 to N, the problem thenbecomes isomorphic to an m-ary search tree. The k-th measurement is usedto determine the k-th m-ary digit in the representation of the channelnumber. Thus, mask 1330 has a vertical layer for each measurement (andhence, sensor) required by the encoding. Each layer subdivides eachregion of the layer above into m sub-regions, each with one of m levelsof transmission.

FIG. 14 is a plot of an exemplary mask 1400 for a 256-channel wavemeterwith binary sensors, in accordance with an embodiment of the presentinvention. White regions are completely transparent, while black regionsare completely opaque. Note that this is an example of the special casefor m=2, where the spectral channel can be determined but not the power.

To work with non-binary sensors requires multiple levels oftransparency. Producing masks with “grayscale” transparency is moredifficult than producing simple binary masks. However, the finitevertical height of individual mask rows can be used to modulate thetransmission to intermediate values, while still using only fullytransmissive or fully opaque coatings.

FIG. 15 is a plot of an exemplary mask 1500 for a 256-channel wavemeterwith 4-level sensors, in accordance with an embodiment of the presentinvention. Again, white regions are completely transmissive and blackregions are completely opaque. The top line of mask 1400 providesinformation about the signal power, and also provides a reference levelfor the remaining rows.

5.2.3 Experimental Results

Another embodiment of the present invention is an exemplary method ofbuilding wavemeter 1300 shown in FIG. 13 for use with mask 1400 of FIG.14 and mask 1500 of FIG. 15. Wavemeter 1300 has a resolutionsubstantially equivalent to one nanometer. Mask 1400 and mask 1500 aredesigned to resolve 256 wavelengths with either four or eightmeasurements.

The resolution and the wavelength range of wavemeter 1300 is determinedby the smallest feature size that can be implemented in either mask andalso by the focal length ‘f’ of the lenses and the grating period ‘Λ’.For a grating with period Λ, the primary grating equation fordiffraction order m is $\begin{matrix}{{{\sin\quad\theta_{i}} + {\sin\quad\theta_{d}}} = \frac{m\quad\lambda}{\Lambda}} & (34)\end{matrix}$where θ_(i) and θ_(d) are the incidence and diffraction anglesrespectively. For a given center wavelength of operation λ_(c), thisequation determines the operating diffraction angle at a particularincidence angle. The grating used has 600 lines per mm and the wavemeteris designed for a center wavelength of 1520 nm. With these values andfor normal incidences, the diffraction angle is 65.8°.

The angular dispersion of the grating is found by differentiating Eq. 34$\begin{matrix}{{\Delta\quad\theta_{d}} = \frac{m\quad{\Delta\lambda}}{\Lambda\quad\cos\quad\theta_{d}}} & (35)\end{matrix}$and the linear dispersion of the grating on the mask is given byΔx=fΔθ_(d)  (36)where f is the focal length of the lens in the 4-f imaging system. Thesedispersion equations determine the feature size of the masks used. A 4-fimaging system uses 50 mm focal length lenses.

An Agilent Tunable laser source 81618A, for example, is used as theinput test source 1305. Laser source 1305 has an operating wavelengthrange of 1466-1586 nm. The range is further restricted to 1470-1570 nm(i.e., 100 channels with a 1 nm even though the masks can resolve up to256 channels). An Indigo Systems AlphaNIR camera, for example, is usedas detector array 1350. Here, the NIR camera is used as only an8-element linear detector measuring only the intensity on each detector.

5.2.4 Binary Mask Results

Another embodiment of the present invention is wavemeter 1300 of FIG. 13including binary mask 1400 of FIG. 14. An Agilent tunable laser source,for example, is used as the test source 1305. The output from the lasersource is tuned to each of the 100 wavelengths in the 1471-1570 nmrange. The fiber coupled output from the source 1305 is collimated and avertical slit of size 25 μm is placed after the collimator. Slit 1310 isimaged on to the mask patterns of mask 1300 using a 4-f imaging systemand through a diffraction grating, for example. Due to grating 1320, theslit that is imaged onto mask 1400 falls at different locations andhence gets spectrally coded differently based on the input wavelength.These patterns uniquely determine the wavelength of input source 1305.

FIG. 16 is a depiction of an exemplary sample image 1600 taken usingexemplary optical wavemeter 1300 of FIG. 13 and the exemplary mask 1400of FIG. 14, in accordance with an embodiment of the present invention.For diagnostic purposes, the output plane of wavemeter 1300 is imagedrather than letting it strike linear detector array 1350. The clearrepresentation of the bit-pattern “10001100” is seen image 1600. Theanalysis of the data is done by simulating the presence of a lineararray in the detection plane. That is, for each simulated photodetector,a pixel region is designated in the image. Further, these regions are ofidentical size and shape, and are regularly spaced along a line.

FIG. 17 is a plot of exemplary expected bit patterns 1700 for each of100 wavelengths using exemplary optical wavemeter 1300 of FIG. 13 andexemplary mask 1400 of FIG. 14, in accordance with an embodiment of thepresent invention. Again, the 100 nm (1471-1570 nm) scan range is due tothe limited tunability of source 1305.

FIG. 18 is a plot of exemplary measured bit patterns 1800 for each of100 wavelengths using the exemplary optical wavemeter 1300 of FIG. 13and the exemplary mask 1400 of FIG. 14, in accordance with an embodimentof the present invention. The 100 wavelength channels are scanned andthe bit patterns obtained on each of the 8 detector element regions onthe NIR camera are measured. 82 unique patterns are observed from themeasurements.

FIG. 19 is a plot of exemplary measured errors 1900 in bit patterns foreach of 100 wavelengths using exemplary optical wavemeter 1300 of FIG.13 and exemplary mask 1400 of FIG. 14, in accordance with an embodimentof the present invention. Measured errors 1900 are the absolute valuesof the error patterns. Errors are chiefly observed in the 7th and 8thbits because of the small feature sizes involved in mask 1400. Also,some errors are caused by the well-known curvature of a slit imagedthrough a grating. This slight curvature is seen in image 1600 of FIG.16 too. As a result, the light (especially in the high resolution areasof the mask) is masked as if it is in a neighboring spectral channel.These errors can be corrected in modifying mask 1400 so that it has aslight arc to match the expected curvature of the slit image.

5.2.4 Base 4 Mask Results

Another embodiment of the present invention is wavemeter 1300 of FIG. 13including 4-level mask 1500 of FIG. 15. Mask 1500 is also designed for256 wavelength channels. Again a laser source 1305 restricts testingonly 100 wavelength channels (the 1471-1570 nm range).

FIG. 20 is a plot of exemplary quaternary digit patterns 2000 for eachof 100 wavelengths using exemplary optical wavemeter 1300 of FIG. 13 andexemplary mask 1500 of FIG. 15, in accordance with an embodiment of thepresent invention. Mask 1500 is a base-4 encoded mask. The values areindicated by the colorbar. The calibration row is suppressed in FIG. 20.

FIG. 21 is a plot of exemplary measured quaternary digit patterns 2100for each of 100 wavelengths using exemplary optical wavemeter 1300 ofFIG. 13 and exemplary mask 1500 of FIG. 15, in accordance with anembodiment of the present invention. The measured intensity at eachsensor is assigned to either 3,2,1 or 0 based on a threshold relative tothe total power as determined by an additional 5th sensor. The valuesare indicated by the colorbar. The calibration row is suppressed in FIG.20.

FIG. 22 is a plot of exemplary measured errors 2200 in quaternary digitpatterns for each of 100 wavelengths using exemplary optical wavemeter1300 of FIG. 13 and exemplary mask 1500 of FIG. 15, in accordance withan embodiment of the present invention. Measured errors 2200 are theabsolute values of the error patterns. Again measured errors 2200 aredue to the curvature, which could be corrected by compensating for it inthe encoding mask.

5.3 Compressive Sampling of Temporal Signals

Another embodiment of the present invention is a method for compressivesampling of a temporal signal. A temporal signal includes but is notlimited to a communications signal, radio wave, or video image pixel.For a temporal signal, measurements are taken asm[n]=s(t),h _(n)(t)

.  (37)

In this case, h_(n)(t) is, for example, implemented electrically. Sincean electrical circuit can be designed to multiply the input signal byany function h_(n)(t), arbitrary compressive sampling is implemented fortemporal signals using an electrical circuit.

FIG. 23 is a flowchart showing an exemplary method 2300 for temporallycompressively sampling a signal, in accordance with an embodiment of thepresent invention.

In step 2310 of method 2300, a plurality of analog to digital convertersis assembled to sample the signal. Each analog to digital converter ofthe plurality of analog to digital converters is configured to samplethe signal at a time step determined by a temporal sampling function.

In step 2320, the signal is sampled over a period of time using theplurality of analog to digital converters. Each analog to digitalconverter of the plurality of analog to digital converters produces ameasurement resulting in a number of measurements for the period oftime.

In step 2330, the number of estimated signal values is calculated fromthe number of measurements and the temporal sampling function. Thetemporal sampling function is selected so that the number ofmeasurements is less than the number of estimated signal values.

6 Compressive Imaging

Another embodiment of the present invention is a method for compressiveimaging. For an image, compressive sampling consists of measuringprojections of the formm[n _(x) ,n _(y) ]=f(x,y),h _(n) _(x,) _(n) _(y) (x,y)

  (38)

The challenge in an imaging system is to obtain non-local response inh_(n) _(x,) _(n) _(y) (x,y). Method for programming h_(n) _(x,) _(n)_(y) (x,y) include but are not limited to setting the shape andinterconnectivity of electronic pixels on a focal plane and/or bymasking pixels, replicating the input image multiple times f(x,y) overredundant apertures and by sampling these resulting images withdifferent focal plane and/or mask structures, and using optical fan-outthrough birefringent structures, micro-optical-mechanical devices ordiffraction gratings. Through a combination of these three mechanismsand linear operations on the output samples, it is possible to implementarbitrary image sampling functions.

6.1 Background

Another embodiment of the present invention is a parallel bank ofimaging systems that is used to synthesize a single global image. Theprimary motivating factor for this embodiment is a reduction in the formfactor of the imaging system. In particular, a system consisting of Nsub-apertures rather than a single coherent aperture is approximately Ntimes thinner than a single lens system. Using compressive sampling, itis possible to also reduce the transverse area of such systems incomparison to non-encoded imagers.

More generally, focal plane coding and compressive sampling enablesuperior data acquisition and handling in imaging systems. Physicallayer data compression reduces the data transfer requirements off of thefocal plane without reducing image resolution. This may be applied inachieving imaging systems with diffraction limited image resolutionwithout electronic sampling at the diffraction limit. Focal plane codingand compressive sampling imaging also enables superior color management,including multispectral and hyperspectral imaging at reasonable datarates. Focal plane coding and compressive sampling can also be used indeveloping polarization cameras. Focal plane coding and compressivesampling imaging strategies can also be applied in temporal focal planesampling to improve data rates and light sensitivity in video imagingsystems.

The focal plane is the interface between the optical field and digitizedimage data. The theory of analog-digital conversion across suchinterfaces began half a century ago with the Whitakker-Shannon samplingtheorem. From the perspective of optical imaging, Whittaker-Shannonindicates that image resolution is proportional to the focal planesampling rate. A quarter century ago, generalized sampling theory showedthat it was possible to reconstruct band-limited signals by sub-Nyquistsampling of parallel multi-channel linear transformations. Manyinvestigators have applied sampling theory to digital imagesuper-resolution.

Focal plane coding and compressive sampling imaging is based on twoimprovements to digital super-resolution techniques. The firstimprovement is sampling function optimization. While the importance ofthe sampling functions to reconstruction fidelity has been clear fromthe origins of generalized sampling, only translation, rotation, andscale have been realized for image super-resolution.

Another embodiment of the present invention is a method of usingmultiplex focal plane encoding to optimize linear transformations forimage reconstruction. Focal plane coding is chosen over aperture orwavefront coding, because previous experience with multiplex imagingsystems has shown that for conventional imaging tasks the highestentropy measurement states are obtained from compact impulse responsetransformations of the focal plane. An emphasis on adaptation of imagingsystem design to the statistics of specific imaging tasks is outlined inthe second improvement.

The second improvement is compressive sampling. Fourier-space samplingtheorems assume that information is uniformly distributed across theimage signal space. The success of image compression algorithms over thepast decade has emphatically demonstrated that this assumption isincorrect. Recognition that image data is spatially and, in the case ofvideo, temporally correlated on multiple scales enables dramaticcompression, particularly with the development of multiscale fractal andwavelet compression bases. Using delocalized multi-scale pixelresponses, focal plane coding enables direct integration of compressivecoding in the image acquisition layer, thereby enabling dramaticimprovements in the data and computational efficiency of digitalsuper-resolution. Another embodiment of the present invention is anoptical sensor employing compressive sampling that tracks a source overN pixels using only log(N) electronic samples.

6.2 Imaging Model

In a model conventional digital imaging system, the transformationbetween the optical image f(x,y) and the focal plane data g_(ij) isg _(ab) =∫∫∫∫p(x−aΔ,y−Δ)h(x−x′,y−y′)f(x′,y′)dxdydx′dy′  (39)where p(x,y) is the pixel sampling distribution on the focal plane, Δ isthe pixel spacing and h(x,y) is the imaging system point spreadfunction.

Another embodiment of the present invention is a method for measuringmultiple demagnified copies of f(x,y) and integrating computationallythese copies into the original full resolution image. Multiple imagecopies are gathered by measuring the same object over multiple imagingsystems, typically using a lenslet array. The final integrated signal isreferred to as “the image” and the copies gathered during themeasurement process are referred to as “subimages.” A conventionalimaging system has an aperture for collecting the field. In thisembodiment the aperture is segmented into multiple “subapertures” foreach subimaging system.

Measurements under this embodiment take the formg _(ab) ^(k) =∫∫∫∫p _(abk)(x,y)h_(k)(x,x′,y,y′)f(mx′,my′)dxdydx′dy′  (40)where k refers to the subimage gathered over the k_(th) subaperture.p_(abk)(x,y) is the measurement sampling function for the ab^(th) pixelin the k^(th) subimage. h_(k)(x,x′,y,y′) is the point spread functionmapping from the density at (x′,y′) in the image to the opticalintensity at (x,y) on the focal plane in the k^(th) subimage. Incontrast with the conventional imaging system, this embodiment allowsfor the possibility that the pixel sampling function and the pointspread functions are shift variant. Shift variance is important increating non-redundant codes between subimages.

The subimages must be non-redundant to enable well-conditioned fullresolution reconstruction of the image f(x,y). Non-redundant codingbetween sub-images is, for example, achieved by the followingmechanisms.

One mechanism involves employing image shifts. Simple shifts in theoptical image relative to the focal plane from one subaperture to thenext make the system weakly non-redundant. This approach is used inTOMBO imaging and similar systems. The TOMBO system demonstratesnondegenerate coding of multiple subimages but does not demonstratecompressive sampling of signal inference. Variations in the registrationof the pixel sampling function relative to the subimage from subapertureto subaperture is the core of this mechanism.

Another mechanism involves using birefringent shift filters.Birefringent shift filters allow controlled and self-registering imageshifts and also allow multiplexed image acquisition and reconstruction.Modulation of the point spread function from subimage to subimage is atthe core of this mechanism.

Another mechanism involves using sampling masks. Sampling masks areplaced in the focal plane to vary p(x,y) from one sub-image to the next.Variations in the spatial structure of the pixel sampling function fromsubimage to subimage is at the core of this mechanism.

Another mechanism involves using microlenslet and microprism arrays.Microlenslet and microprism arrays are be used to remap the opticalfocal plane from subimage to subimage. This mechanism modulates thepoint spread function in a shift-variant manner.

Another mechanism involves using spectral and spatial filters. Spectralfilters, such as gratings and Fabry-Perot and Fizeau Filters can remapthe spatio-spectral pattern of the image from one sub-aperture to thenext. This mechanism modulates the structure of the point spreadfunction, typically in a chromatically varying fashion.

Another mechanism involves the physical structure of the focal plane andfocal plane read dynamics. The physical structure of the focal plane andfocal plane read dynamics may vary from one subimage and from one pixelto the next. This mechanism is similar to foveation. This mechanismmodulates the structure of the pixel sampling function more profoundlythan can be done using simple pixel masks.

Another mechanism involves using micro-optical mechanical devices.Micro-optical mechanical devices, such as deformable mirror arrays,digital light projectors, and similar dynamic elements are used toencode the optical impulse response for compressive sampling. Liquidcrystal or electro-optic spatial light modulators might are similarlyapplied to effectively combine the functions of coding masks andmicro-optical elements. Dynamic spatial light modulators are used toadapt sensor encodings to improve compressive sampling.

6.3 Coding Mechanisms

Another embodiment of the present invention is a focal plane coding andcompressive sampling imaging system that measures linear transformationsof an image f(x,y) such that the full image is digitally reconstructed.This embodiment uses optical and electro-optical prefilters to encodethe linear transformations to enable and improve the fidelity of imagereconstruction.

While in some cases this embodiment is used in single aperture systems,in most cases this embodiment involves encoding over multiplesubapertures. In this case, the imaging system is regarded as a filterbank and/or sampling system array.

6.3.1 Pixel Shift Coding

Another embodiment of the present invention is a method for pixel shiftcoding. This method can be used individually to encode images forreconstruction or can be used in combination with other methods. Pixelshift coding consists of measurement under different pixel shifts fromone subaperture to the next. The measurement model for this case isg _(ab) ^(k) ^(x) ^(k) ^(y) ∫∫∫∫p(x−aΔ−k _(x) δ,y−bΔ−k_(y)δ)h(x−x′,y−y′)f(Mx′,My′)dxdydx′dy′  (41)where δ is a subimage to subimage differential shift. While it isassumed here that there is a constant differential shift, thedifferential pixel shift, in practice, can vary nonlinearly and imagerotations can also be considered from one subimage to the next.

Pixel shift coding has been previously considered by many. Some, inparticular, have considered the use of pixel shift coding using multiplesubimage arrays. Registration of the subimages relative to the globalimage coordinates has been a primary challenge for these methods. Othershave focused on interpolation of full images based on self-registrationusing correlations between captured images.

The total pixel shift from subimage to subimage is due to shifts in theregistration of the optical axis of each subaperture relative to theelectronic pixel axis. If, for example, the subapertures are arranged ona Cartesian grid with uniform spacing Δ_(lens) between optical axis inboth x and y, the subimage differential shift isδ=mod(Δ_(lens),Δ)  (42)where mod( ) is the modulus function.

The total pixel shift from subimage to subimage is also due to shifts inthe registration of the optical axis of each subaperture relative to theaxis of the object. The differential shift due to the subaperture-objectshift is $\begin{matrix}{\delta = {\Delta_{lens}\frac{F}{R}}} & (43)\end{matrix}$where F is the focal length of the imaging lenslet and R is the objectrange.

The object range dependent differential shift makes it difficult toimplement uniform coding using pixel shifts and requires theself-registration techniques used by several authors.

The value of the differential pixel shifts critically determines theconditioning of the subimages for full image estimation. Uniformlyspaced shifts spanning the range of the pixel sampling period areparticularly attractive. The simplest approach to understanding imageestimation for shifted pixel imagers relies on expansion of p(x,y) on afinite basis. In typical systems, p(x,y) is a uniform rectangularfunction with finite fill factor. Assuming 100% fill factor it can beassumed that $\begin{matrix}{{p\left( {x,y} \right)} = {\beta^{o}\left( \frac{x}{\Delta} \right)}} & (44)\end{matrix}$where β⁰( ) is the order 0 B-spline. If Δ/δ is an integer, the family offunctions $\begin{matrix}{{\phi_{n}(x)} = {\beta^{o}\left( \frac{x - {n\quad\delta}}{\delta} \right)}} & (45)\end{matrix}$spans the range of the pixel sampling functions specifically$\begin{matrix}{{p\left( {{x - {a\quad\Delta} - {k_{x\quad}\delta}},{y - {b\quad\Delta} - {k_{y}\delta}}} \right)} = {\sum\limits_{n = {{{({a - \frac{1}{2}})}\frac{\Delta}{\delta}} + k_{x} + \frac{1}{2}}}^{n = {{{({a + \frac{1}{2}})}\frac{\Delta}{\delta}} + k_{x} - \frac{1}{2}}}\quad{\sum\limits_{n = {m = {{{({b - \frac{1}{2}})}\frac{\Delta}{\delta}} + k_{y} + \frac{1}{2}}}}^{n = {{{({b + \frac{1}{2}})}\frac{\Delta}{\delta}} + k_{y} - \frac{1}{2}}}\quad{{\phi_{n}(x)}{\phi_{m}(y)}}}}} & (46)\end{matrix}$where it is assumed that Δ/δ is even. Defining a canonicalrepresentation for the discrete image f_(c) asf _(c)(n,m)=∫∫∫∫φ_(n)(x)φ_(m)(y)h(x−x′,y−y′)f(Mx′,My′)dxdydx′dy′  (47)Eqn. (41) is expressed in discrete form asg=Tf_(c)  (48)where T consists of N² ap submatricies, one for each subimage. Each rowof T contains)Δ/δ² ones and all other components are zero. Although T isnot particularly well conditioned, f_(c) is estimated by directalgebraic inversion or other methods.

A critical challenge for pixel shift methods, is that Δ/δ is notgenerally an integer, in which case there is no convenient finite basisfor representation of f_(c). In this case, interpolation methods arerequired.

6.3.2 Birefringent Shift Coding

Another embodiment of the present invention is a method for birefringentshift coding. Birefringent coding elements having previously been usedto produce coded optical fan-out. Each layer of a birefringent filterproduces two orthogonally polarized copies of the incident field. It isassumed that the image field is unpolarized, although coding to estimatepolarization image is also possible. The transformation with shiftcoding induced by a birefringent filter is $\begin{matrix}{g_{ab}^{k} = {\int{\int{\int{\int{{p\left( {{x - {a\quad\Delta}},{y - {b\quad\Delta}}} \right)}{h\left( {{x - x^{\prime}},{y - y^{\prime}}} \right)}{\sum\limits_{n}\quad{{f\left( {{{M\quad x^{\prime}} - {\delta\quad x_{nk}}},{{My}^{\prime} - {\delta\quad y_{nk}}}} \right)}{\mathbb{d}x}{\mathbb{d}y}{\mathbb{d}x^{\prime}}{\mathbb{d}y^{\prime}}}}}}}}}} & (49)\end{matrix}$

This transformation is equivalently expressed in terms of the Fouriertransforms of the pixel sampling function P(u,v), the point spreadfunction (PSF) H(u,v), and the intensity distribution F(u,v) as$\begin{matrix}{g_{ab}^{k} = {\frac{1}{M^{2}}{\int{\int{{P\left( {u,v} \right)}{H\left( {u,v} \right)}{F\left( {\frac{u}{M},\frac{v}{M}} \right)}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad\pi\quad u\quad a\quad\Delta}{{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad\pi\quad v\quad b\quad\Delta}\left\lbrack {\sum\limits_{n}{{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad\pi\quad u\quad\delta\quad x_{nk}}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad\pi\quad v\quad\delta\quad y_{nk}}}} \right\rbrack}{\mathbb{d}u}{\mathbb{d}v}}}}}} & (50)\end{matrix}$

The discrete Fourier transform (DFT) of the measured data is:$\begin{matrix}\begin{matrix}{G_{\alpha\quad\beta}^{k} = {\frac{1}{N_{S}}{\sum\limits_{{ab} = 0}^{N_{s} - 1}{g_{ab}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad\pi\quad\frac{a\quad\alpha}{N_{s}}}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad\pi\quad\frac{b\quad\beta}{N_{s}}}}}}} \\{= {\frac{1}{N_{s}M^{2}}{\int{\int{{P\left( {u,v} \right)}{H\left( {u,} \right.}}}}}} \\{\left. v \right){F\left( {\frac{u}{M},\frac{v}{M}} \right)}{\mathbb{e}}^{{- {\mathbb{i}}}\quad{\pi{({N_{s} - 1})}}{({\frac{\alpha}{N_{s}} + {u\quad\Delta}})}}{\mathbb{e}}^{{- {\mathbb{i}}}\quad{\pi{({N_{s} - 1})}}{({\frac{\beta}{N_{s}} + {v\quad\Delta}})}} \times} \\{\frac{\sin\left( {\pi\left( {\alpha + {N_{s}u\quad\Delta}} \right)} \right)}{\sin\left( {\pi\left( {\frac{\alpha}{N_{s}} + {u\quad\Delta}} \right)} \right)}\frac{\sin\left( {\pi\left( {\beta + {N_{s}v\quad\Delta}} \right)} \right)}{\sin\left( {\pi\left( {\frac{\beta}{N_{s}} + {v\quad\Delta}} \right)} \right)}} \\{\left\lbrack {\sum\limits_{n^{\prime}}{{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\pi\quad u\quad\delta\quad x_{n^{\prime}k}}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad\pi\quad v\quad\delta\quad y_{n^{\prime}k}}}} \right\rbrack{\mathbb{d}u}{\mathbb{d}v}} \\{\approx {\frac{1}{N_{s}M^{2}}{\sum\limits_{n = {- \infty}}^{\infty}{\sum\limits_{m = {- \infty}}^{\infty}{P\left( {{\frac{n}{\Delta} - \frac{\alpha}{\Delta\quad N_{s}}},{\frac{m}{\Delta} -}} \right.}}}}} \\{\left. \frac{\beta}{\Delta\quad N_{s}} \right){H\left( {{\frac{n}{\Delta} - \frac{\alpha}{\Delta\quad N_{s}}},{\frac{m}{\Delta} - \frac{\beta}{\Delta\quad N_{s}}}} \right)} \times} \\{F\left( {{\frac{n}{M\quad\Delta} - \frac{\alpha}{M\quad\Delta\quad N_{s}}},{\frac{m}{M\quad\Delta} - \frac{\beta}{M\quad\Delta\quad N_{s}}}} \right)} \\{\left\lbrack {\sum\limits_{k}{{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad{\pi{({\frac{n}{\Delta} - \frac{\alpha}{\Delta\quad N_{s}}})}}\delta\quad x_{k}}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad{\pi{({\frac{m}{\Delta} - \frac{\beta}{\Delta\quad N_{s}}})}}{\delta y}_{k}}}} \right\rbrack}\end{matrix} & (51)\end{matrix}$

Assuming full fill factor rectangular pixels, P(u,v)=sinc(uΔ)sinc(uΔ).For a diffraction limited imaging H(u,v)=β¹(uλ/2N.A., uλ/2N.A.), whereN.A. is the numerical aperture. β¹ is the first order B-spline.

Assuming that α and β are limited to [0, N_(s)−1], H(u,v) isnonvanishing in the sum of Eqn. (51) for: $\begin{matrix}{{n \leq {\frac{\alpha}{N_{s\quad}} + \frac{2\quad\Delta\quad{N.A.}}{\lambda}}}{m \leq {\frac{\beta}{N_{s}} + \frac{2\Delta\quad{N.A.}}{\lambda}}}} & (52)\end{matrix}$This means that approximately$N_{sub} = {4\left( \frac{\Delta\quad{N.A.}}{\lambda} \right)^{2}}$

independent samples of F(u,v) are multiplexed in each Fourier sample.Ideally, N_(sub) copies of the image are sampled with differentbirefringent shifts to invert and reconstruct.

In the simplest case, it is assumed that just two shifted copies of theimage in each subaperture. In this case, $\begin{matrix}{\left\lbrack {\sum\limits_{k}{{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad{\pi{({\frac{n}{\Delta} - \frac{\alpha}{\Delta\quad N_{s}}})}}\delta\quad x_{k}}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad{\pi{({\frac{n}{\Delta} - \frac{\alpha}{\Delta\quad N_{s}}})}}\delta\quad y_{k}}}} \right\rbrack = {\cos\begin{pmatrix}{{2{\pi\left( {\frac{n}{\Delta} - \frac{\alpha}{\Delta\quad N_{s}}} \right)}\delta\quad x_{k}} +} \\{2{\pi\left( {\frac{m}{\Delta} - \frac{\alpha}{\Delta\quad N_{s}}} \right)}\delta\quad y_{k}}\end{pmatrix}}} & (53)\end{matrix}$The set of measurements over the subapertures of$F\left( {{\frac{n}{M\quad\Delta} - \frac{\alpha}{M\quad\Delta\quad N_{s}}},{\frac{m}{M\quad\Delta} - \frac{\beta}{M\quad\Delta\quad N_{s}}}} \right)$is a DCT of Fourier coefficients over the range of n and m. Ideallyinversion of this DCT requires $\begin{matrix}{{\delta_{x} = \frac{n^{\prime}\Delta}{\sqrt{N_{sub}}}}{\delta_{y} = \frac{m^{\prime}\Delta}{\sqrt{N_{sub}}}}} & (54)\end{matrix}$where (n′,m′) is the subaperture index.

The pixel shift coding can be analyzed by similar Fourier methods.However, birefringent coding offers several improvements over pixelshift coding. In particular, birefringent shifts are coded by settingthe thickness and orientation of a birefringent plate. Birefringentfan-out can be programmed at the tolerances required for thisapplication. Exact coding of shifts has not been demonstrated to thesame tolerance.

Also, birefringent shifts are self-registering. The challenge of shiftcoding is that the shift is object range dependent and alignmentdependent. Subimage shifts in a birefringently coded system cause aglobal phase shift on the subimage Fouier transform. Given knownconstraints on the image (i.e., it is real and nonnegative), previousstudies have shown that the image can be recovered from the magnitude ofits Fourier transform. Thus, one algorithm can eliminate shift varianceby considering |G_(αβ) ^(k)|. This system still contains the cosinetransform of the modulus of F(u,v). In the birefringent coding case, theimage is coded against its own reference.

6.3.3 Transmission Mask Coding

Another embodiment of the present invention is a method for birefringentshift coding. The transformation for coding with a focal plane mask is:g _(ij) =∫∫t(x,y)p(x−iΔ, y−jΔ)h(x−x, y−y′)f(mx′, my′)dxdydx′dy′  (55)The idea of coding with a mask is to select t(x,y) such that the lineartransformation is well-conditioned for conventional or compressivelysampled reconstruction.6.3.3 Microoptic, Diffractive, and Interferometric Coding

While birefringent coding offers the simplest and least dispersiveapproach to optical prefiltering focal plane coding and compressivesampling encoding, micro-optical elements, such as lenslets, diffractiveelements and prisms, can be used to reshape the impulse response in acoded manner for compressive or non-compressive imaging. Interferometersin the optical path, including Fabry-Perot, Fizeau and Mach-Zenderdevices, can also code the impulse response.

Another embodiment of the present invention is a method for opticalprefilter encoding. The basic form of the transformation for opticalprefilter encoding isg _(ij) ^(k) =∫∫p(x−iΔ, y−Δ)h _(k)(x,x′,y,y′,λ)f(mx′,my′λ)dxdydx′dy′dλ  (56)with the basic point being that the point spread function is not shiftinvariant and that it is non-redundantly coded in each subaperture. Thepossibility that the point spread function is dispersive (wavelengthdependent) is also include.

Wavelength dispersion is both a challenge and an opportunity. In casesfor which one is more interested in the achromatic image, wavelengthdispersion can decrease signal to noise by binning light from similarspatial location in different detection channels. In cases where one isinterested in multispectral or hyperspectral images, wavelengthdispersion offers an opportunity to create a coded cross-spectralmeasurement system.

Temporal variations and adaptation in the optical prefilter is anotherpossibility with micro-optical signal encoding. Using a spatial lightmodulator, such as a micromirror array, a liquid crystal or acousticmodulator, it is possible to vary h_(k)(x,x′,y,y′, λ)f(mx′,my′, λ) as afunction of time. This can be used to implement time domain compressivesampling or to match the PSF to the image signal.

6.4 Focal Plane Coding

Another embodiment of the present invention is a method for focal planecoding. A focal plane is an array of electronic detectors. In typicalapplications the array has a grid of rectangular pixels. Rectangularpixels are particularly appropriate for, for example, charge coupleddevice based detectors (CCDs). Other electronic detectors include activepixel CMOS VLSI detector arrays. In infrared systems, CMOS or CCDsystems are used for the backend or front-end detector array.

Using CMOS detectors, microbolometer, or other detector arrays it ispossible to shape both the spatial distribution of pixel samplingfunctions and the temporal read-out circuit for compressive sampling.The focal plane modulated image to measurement mapping isg _(ij) ^(k)(t)=∫∫p _(k)(x−iΔ,y−jΔ,t)h(x−x′,y−y′)f(mx′,my′)dxdydx′dy′  (57)The possibility of temporal variations is included in the pixel samplingfunction. This constitutes temporal variation in the read-out rate andread-out filter from different pixels rather than temporal changes inthe spatial structure.

Generalized filtering in the temporal read-out circuit enablescompressive sampling in time as well as space. From a spatialperspective, the lay-out of the pixel sampling structure achieves thesame advantages as a transmission mask without masking. In addition,however, electronic interconnections and space-time filtering enablefocal planes to implement non-convex, disjoint and multiscale pixelsampling functions that can be difficult to achieve with simple masksampling.

FIG. 24 is a flowchart showing an exemplary method 2400 for spatiallyfocal plane coding an optical signal using an imaging system, inaccordance with an embodiment of the present invention.

In step 2410 of method 2400, the imaging system is created from aplurality of subimaging systems. Each subimaging system includes anarray of electronic detectors. Each subimaging system is dispersedspatially with respect to a source of the optical signal according to asampling function.

In step 2420, the optical signal is collected at each array ofelectronic detectors of the plurality of subimaging systems at a singletime step.

In step 2430, the optical signal is transformed into a subimage at eachsubimaging system. The subimage includes at least one measurement fromthe array of electronic detectors of the subimaging system.

In step 2440, an image of the optical signal is calculated from thesampling function and each subimage, spatial location, pixel samplingfunction, and point spread function of each subimaging system of theplurality of subimaging systems. The sampling function is selected sothat the number of measurements for the plurality of subimages is lessthan the number of estimated optical signal values that make up theimage.

FIG. 25 is a flowchart showing an exemplary method 2500 for temporallyfocal plane coding an optical signal using an imaging system, inaccordance with an embodiment of the present invention.

In step 2510 of method 2500, the imaging system is created from aplurality of subimaging systems. Each subimaging system includes anarray of electronic detectors.

In step 2520, the optical signal is collected at each array ofelectronic detectors of the plurality of subimaging systems at adifferent time step according to a temporal sampling function.

In step 2530, the optical signal is transformed into a subimage at eachsubimaging system of the plurality of subimaging systems. The subimageincluded at least one measurement from an array of electronic detectorsof each subimaging system.

In step 2540, an image of the optical signal is calculated from thetemporal sampling function and each subimage, time step, pixel samplingfunction, and point spread function of each subimaging system of theplurality of subimaging systems. The temporal sampling function isselected so that the number of measurements for the plurality ofsubimages is less than the number of estimated optical signal valuescomprising the image.

FIG. 26 is a flowchart showing an exemplary method 2600 for spatiallyand temporally focal plane coding an optical signal using an imagingsystem, in accordance with an embodiment of the present invention.

In step 2610 of method 2600, the imaging system is created from aplurality of subimaging systems. Each subimaging system includes anarray of electronic detectors. Each subimaging system is dispersedspatially with respect to a source of the optical signal according to aspatial sampling function.

In step 2620, the optical signal is collected at each array ofelectronic detectors of the plurality of subimaging systems at adifferent time step according to a temporal sampling function.

In step 2630, the optical signal is transformed into a subimage at eachsubimaging system of the plurality of subimaging systems. The subimageincludes at least one measurement from an array of electronic detectorsof the each subimaging system.

In step 2640, an image of the optical signal is calculated from thespatial sampling function, the temporal sampling function and eachsubimage, spatial location, time step, pixel sampling function, andpoint spread function of each subimaging system of the plurality ofsubimaging systems. The spatial sampling function and the temporalsampling function are selected so that the number of measurements forthe plurality of subimages is less than the number of estimated opticalsignal values that make up the image.

7 Dynamic and Adaptive Coding Systems

It is useful in some cases to apply dynamic, rather than static,sampling strategies in compressive sampling systems. In such cases,micro-electro-optical mechanical devices or other dynamic spatial lightmodulators are used to adapt sampling filters and codes.

In accordance with an embodiment of the present invention, instructionsconfigured to be executed by a processor to perform a method are storedon a computer-readable medium. The computer-readable medium can be adevice that stores digital information. For example, a computer-readablemedium includes a compact disc read-only memory (CD-ROM) as is known inthe art for storing software. The computer-readable medium is accessedby a processor suitable for executing instructions configured to beexecuted. The terms “instructions configured to be executed” and“instructions to be executed” are meant to encompass any instructionsthat are ready to be executed in their present form (e.g., machine code)by a processor, or require further manipulation (e.g., compilation,decryption, or provided with an access code, etc.) to be ready to beexecuted by a processor.

Systems and methods in accordance with an embodiment of the presentinvention disclosed herein can advantageously reduce sensor cost, reducesensor data load, reduce processing load, increase spatial samplingcapacity and improve spatial abstraction.

For spectral bands at the edge of the visible spectrum, such asinfrared, millimeter wave and ultraviolet ranges detector cost is asignificant component of overall system cost. Compressive samplingenables lower system cost by reducing the detector count necessary for agiven range and resolution specification.

In highly parallel sensor systems, such as imaging systems, read-outdata loads may be overwhelming. Compression at the point of sampling mayreduce these data loads and enable lower power, lower bandwidth sensordesign.

Many sensor systems implement data compression in software immediatelyafter sensing. Direct compression in sampling reduces or eliminates postsensor processing for data compression and analysis.

Sensor systems are limited by physical geometric constraints.Hyperspectral imagers, for example, require simultaneous sampling ofmyriad spatial and spectral channels. Efficient packing of samplingpoints requires mechanical scanning or large system apertures.Compressive sampling enables spatial and spectral sampling to be foldedinto smaller system volumes and shorter scanning times.

Sensor systems are based on naive isomorphisms between object andmeasurement spaces. The discrete abstractions of compressive samplingenable generalized models for the object space, includingmultidimensional object spaces and differences between thedimensionality of the measurement and object space.

8 Multiplex and Reflective Compressing Sampling

FIG. 27 is a flowchart showing an exemplary method 2700 forcompressively sampling an optical signal using multiplex sampling, inaccordance with an embodiment of the present invention.

In step 2710 of method 2700, group measurement of multiple opticalsignal components on single detectors is performed such that less thanone measurement per signal component is taken.

In step 2720, probabilistic inference is used to produce a decompressedoptical signal from the group measurement on a digital computer.

FIG. 28 is a flowchart showing an exemplary method 2800 forcompressively sampling an optical signal using a reflective mask, inaccordance with an embodiment of the present invention.

In step 2810 of method 2800, an reflective mask with a plurality ofreflective elements and a plurality of non-reflective elements iscreated. The location of the plurality of reflective elements and theplurality of non-reflective elements is determined by a reflectionfunction. The reflective mask includes but is not limited to areflection mask or a coded aperture.

In step 2820, the spectrum of the optical signal is dispersed across thereflective mask.

In step 2830, signals reflected by the plurality of reflective elementsare detected in a single time step at each sensor of a plurality ofsensors dispersed spatially with respect to the reflective mask. Eachsensor of the plurality of sensors produces a measurement resulting in anumber of measurements for the single time step.

In step 2840, a number of estimated optical signal values is calculatedfrom the number of measurements and the reflection function. Thereflection function is selected so that the number of measurements isless than the number of estimated optical signal values.

In another embodiment of the present invention, the number of estimatedoptical signal values is calculated by multiplying the number ofmeasurements by a pseudo-inverse of the reflection function. In anotherembodiment of the present invention, the number of estimated opticalsignal values is calculated by a linear constrained optimizationprocedure. In another embodiment of the present invention, the number ofestimated optical signal values is calculated by a non-linearconstrained optimization procedure. In another embodiment of the presentinvention, the reflective mask is adaptively or dynamically variable. Inanother embodiment of the present invention, the adaptively ordynamically variable reflective mask is based on electro-optic effects,liquid crystals, or micro mechanical devices. In another embodiment ofthe present invention, the reflective mask implements positive andrespective negative elements of an encoding.

FIG. 29 is a flowchart showing an exemplary method 2900 forcompressively sampling an optical signal using an optical component toencode multiplex measurements, in accordance with an embodiment of thepresent invention.

In step 2910 of method 2900, a mapping is created from the opticalsignal to a detector array by spatial and/or spectral dispersion.

In step 2920, signals transmitted are detected by a plurality oftransmissive elements of the optical component at each sensor of aplurality of sensors of the detector array dispersed spatially withrespect to the optical component. Each sensor of the plurality ofsensors produces a measurement resulting in a number of measurements.

In step 2930, a number of estimated optical signal values is calculatedfrom the number of measurements and a transmission function, wherein thetransmission function is selected so that the number of measurements isless than the number of estimated optical signal values.

Another embodiment of the present invention is a spectrometer usingcompressive sampling. The spectrometer includes a plurality of opticalcomponents and a digital computer. The plurality of optical componentsmeasure multiple linear projections of a spectral signal. The digitalcomputer performs decompressive inference on the multiple linearprojections to produce a decompressed signal.

The foregoing disclosure of the preferred embodiments of the presentinvention has been presented for purposes of illustration anddescription. It is not intended to be exhaustive or to limit theinvention to the precise forms disclosed. Many variations andmodifications of the embodiments described herein will be apparent toone of ordinary skill in the art in light of the above disclosure. Thescope of the invention is to be defined only by the claims appendedhereto, and by their equivalents.

Further, in describing representative embodiments of the presentinvention, the specification may have presented the method and/orprocess of the present invention as a particular sequence of steps.However, to the extent that the method or process does not rely on theparticular order of steps set forth herein, the method or process shouldnot be limited to the particular sequence of steps described. As one ofordinary skill in the art would appreciate, other sequences of steps maybe possible. Therefore, the particular order of the steps set forth inthe specification should not be construed as limitations on the claims.In addition, the claims directed to the method and/or process of thepresent invention should not be limited to the performance of theirsteps in the order written, and one skilled in the art can readilyappreciate that the sequences may be varied and still remain within thespirit and scope of the present invention.

1. A method for compressively sampling an optical signal, comprising:creating an optical component with a plurality of transmissive elementsand a plurality of opaque elements, wherein the location of theplurality of transmissive elements and the plurality of opaque elementsis determined by a transmission function; dispersing a spectrum of theoptical signal across the optical component; detecting signalstransmitted by the plurality of transmissive elements in a single timestep at each sensor of a plurality of sensors dispersed spatially withrespect to the optical component, wherein the each sensor of a pluralityof sensors produces a measurement resulting in a number of measurementsfor the single time step; and calculating a number of estimated opticalsignal values from the number of measurements and the transmissionfunction, wherein the transmission function is selected so that thenumber of measurements is less than the number of estimated opticalsignal values.
 2. The method of claim 1, wherein the optical componentcomprises one of a transmission mask and a coded aperture.
 3. The methodof claim 1, wherein the number of estimated optical signal values iscalculated by multiplying the number of measurements by a pseudo-inverseof the transmission function.
 4. The method of claim 1, wherein thenumber of estimated optical signal values is calculated by a linearconstrained optimization procedure.
 5. The method of claim 1, whereinthe number of estimated optical signal values is calculated by anon-linear constrained optimization procedure.
 6. The method of claim 1,wherein the dispersing a spectrum of the optical signal across theoptical component is performed using one of a hologram and a diffractiongrating.
 7. The method of claim 1, wherein the optical component isadaptively variable.
 8. The method of claim 1, wherein the opticalcomponent is dynamically variable.
 9. The method of claim 1, wherein theoptical component comprises one or more of an electro-optic effect, aliquid crystal, and a micro mechanical device.
 10. A spectrometer usingcompressive sampling, comprising: a transmission mask with a pluralityof transmissive elements and a plurality of opaque elements, wherein thelocation of the plurality of transmissive elements and the plurality ofopaque elements is determined by a transmission function; a diffractiongrating, wherein the diffraction grating disperses a spectrum of anoptical signal across the transmission mask; a plurality of sensorsdispersed spatially with respect to the transmission mask, wherein eachsensor of the plurality of sensors detects signals transmitted by theplurality of transmissive elements in a single time step, and whereinthe each sensor of a plurality of sensors produces a measurementresulting in a number of measurements for the single time step; and aprocessor, wherein the processor calculates a number of estimatedoptical signal values from the number of measurements and thetransmission function, and wherein the transmission function is selectedso that the number of measurements is less than the number of estimatedoptical signal values.
 11. The spectrometer of claim 10, wherein thetransmission mask comprises one of a photonic crystal filter, aholographic filter, a holographic filter array, a multiplex gratinghologram, a multipass interferometer, a thin film filter, an array ofthin film filters, and an array of photonic crystal filters.
 12. Amethod for compressively sampling an optical signal, comprising:creating a reflective mask with a plurality of reflective elements and aplurality of non-reflective elements, wherein the location of theplurality of reflective elements and the plurality of non-reflectiveelements is determined by a reflection function; dispersing a spectrumof the optical signal across the reflective mask; detecting signalsreflected by the plurality of reflective elements in a single time stepat each sensor of a plurality of sensors dispersed spatially withrespect to the reflective mask, wherein the each sensor of a pluralityof sensors produces a measurement resulting in a number of measurementsfor the single time step; and calculating a number of estimated opticalsignal values from the number of measurements and the reflectionfunction, wherein the reflection function is selected so that the numberof measurements is less than the number of estimated optical signalvalues.
 13. The method of claim 12, wherein the reflective mask isadaptively variable.
 14. The method of claim 12, wherein the reflectivemask is dynamically variable.
 15. The method of claim 12, wherein thereflective mask comprises one or more of an electro-optic effect, aliquid crystal, and a micro mechanical device.